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EXECUTIVE  SUMMARY 


Recently  the  Navy  has  felt  the  need  for  an  improved  common  tactical  picture  that  incorporates 
information  from  all  sources.  These  include  Navy,  Joint,  Allied,  and  Coalition  sources  as  well  as 
sensor  information,  passive  and  active,  from  air,  surface,  submarine,  and  space  platforms.  In  most 
cases,  the  types  of  information  are  significantly  complex  and  nonlinear.  As  the  Navy  moves  into 
operating  in  more  littoral  regions,  the  sensor  information  will  take  on  an  even  more  complex 
character.  It  has  quickly  become  apparent  that  the  present  methods  and  approaches  to  data  fusion, 
which  are  based  on  Kalman  filtering  and  its  extensions,  are  inadequate  to  cope  with  this  problem. 


The  Nodestar  engine,  described  in  this  report,  is  extremely  powerful  and  represents  an  advance  in 
the  state  of  the  art  for  data  fusion  and  decision  support.  Additionally,  extensions  to  Nodestar  are 
straightforward.  For  any  observation  or  response  representing  a  new  type  of  information  or  sensor, 
one  must  simply  construct  a  likelihood  function  that  computes  the  probability  of  receiving  that 
response  as  a  function  of  target  state.  This  usually  requires  two  things:  first,  a  good  physical  model 
of  the  sensor  or  response  process  and,  second,  that  the  target  state  space  contain  the  right  information 
for  the  calculation,  e.g.,  velocity,  if  one  wishes  to  account  for  Doppler  in  a  frequency  reading.  If  the 
state  space  contains  the  necessary  information,  then  one  has  only  to  add  a  new  module  for  calculating 
the  likelihood  function.  Once  this  is  done,  new  information  can  be  incorporated  via  the  information 
update  function. 

We  describe  the  implementation  of  Nodestar  in  the  Spotlight  decision  support  system.  The 
Spotlight  system  demonstrated  an  innovative  computer-based  approach  to  managing  and  developing 
resource  optimization  plans  for  a  deep  water  ASW  surveillance  post.  Extensions  to  surface  ship 
tracking  and  active  sonar  tracking  are  described  as  well. 

The  bulk  of  this  work  was  performed  under  Naval  Research  Laboratory  (NRL)  contract  N00014- 
91-C-2140.  Development  of  the  Nodestar  nonlinear  discrete  tracker  was  begun  at  NRL  in  1988  as 
part  of  the  Tactical  ASW  Battle  Management  System  (TABS)  under  the  Command  Decision  Support 
6.2  program.  Nodestar  reached  its  culmination  in  the  multiple-target  tracker  developed  for  the 
Spotlight  Advanced  Technology  Demonstration  (ATD).  An  extension  of  the  algorithm,  also 
documented  herein,  was  created  for  multistatic  active  tracking  with  the  Shallow  Water  Active  Tracking 
System  (SWATS).  SWATS  work  has  been  performed  as  part  of  a  Phase  I  SBIR  from  NAVSEA  and 
continues  under  Phase  II  funding.  Appendix  A  summarizes  the  technical  background  for  modeling 
the  spherical  motion.  A  summary  of  the  performance  comparison  of  the  Nodestar  nonlinear  tracker 
(Version  1991)  with  MTST  (Maneuvering  Target  Statistical  Tracker)  linear  tracker  appears  in 
Appendix  B. 

The  Spotlight  ATD  included  the  development  of  a  multiple-target,  nonlinear  correlator-tracker 
(Nodestar)  and  a  Resource  Optimization  (RO)  system  for  the  Integrated  Underwater  Surveillance 
System  (lUSS).  The  correlator-tracker  and  the  RO  system  in  Spotlight  are  integrated  through  an 
innovative  Operator-Machine  Interface  (OMI).  The  RO  system  is  documented  in  a  separate  report. 
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Nodcstar  was  demonstrated  as  part  of  the  Spotlight  ATD  test  in  December  1993  at  the  Dam  Neck 
Naval  Ocean  Processing  Facility. 

The  Spotlight  version  of  Nodestar  processes  line-of-bearing  information  (including  mirror 
bearings)  and  posits  (ellipses)  while  incorporating  performance  prediction  information  from  System 
for  the  Prediction  of  Acoustic  Responses  (SPARSE  It  processes  negative  (detection/no-detection) 
information  and  correctly  accounts  for  the  presence  of  land.  The  state  space  is  six-dimensional  with 
the  following  components:  position  (3D),  velocity  (2D),  and  target  class  (ID).  Nodestar  uses 
information  about  target  class,  expressed  in  terms  of  probability  distributions,  to  help  associate  and 
fuse  frequency  information  contained  in  detection  reports.  This  process  requires  the  use  of  a  threat 
data  base  that  identifies  the  set  of  frequencies  radiated  by  each  target  class.  This  information  in  turn 
is  used  to  update  the  distribution  on  the  target's  class.  A  similar  process  is  followed  to  process 
information  about  electronic  intelligence  (FLINT)  parameters. 


In  comparison  to  the  restrictive  assumptions  inherent  in  most  current  Kalman-filter  based  tracker- 
correlators,  the  ones  made  for  Nodestar  allow  for  much  more  freedom  and  fidelity  in  modeling  target 
motion  and  sensor  response  and  as  such  represent  an  advance  in  the  state  of  the  art  in  data  fusion. 
The  unique  features  of  Nodestar  were  exploited  throughout  the  Spotlight  system  to  perform  more 
efficiently  the  functions  of  forecasting,  “what  if?”  analysis,  time-late  data  inclusion,  data  correction, 
and  the  aforementioned  Resource  Optimization  module.  As  implemented  in  Spotlight,  Nodestar  also 
represents  an  advance  in  the  state  of  the  art  for  decision  support. 


E-2 


TECHNICAL  DOCUMENTATION  OF  NODESTAR 


1.  INTRODUCTION 

In  this  report,  we  present  a  framework  in  which  significantly  complex  and  nonlinear  sensor  data 
can  be  correctly  and  effectively  addressed.  We  describe  a  general  framework  for  the  process  of 
multiple  target,  multiple  sensor  data  fusion.  With  that  framework  in  place,  we  identify  those  areas  in 
which  the  methodology  is  mature  and  ripe  for  implementation  and  those  areas  that  require  further 
development.  One  area  in  which  the  methodology  is  well  advanced  is  nonlinear  tracking.  For  that 
area  we  describe  a  basic  engine.  Nodestar,  which  has  been  developed  to  perform  nonlinear,  multiple- 
target  tracking.  Then  we  describe  the  version  of  Nodestar  that  has  been  developed  for  the  Spotlight 
Advanced  Technology  Demonstration.  We  finish  with  a  discussion  of  extensions  to  this  version  of 
Nodestar. 

Using  the  Nodestar  engine  one  can  incorporate  new  types  of  information  by  the  process  of 
calculating  likelihood  functions.  Once  the  likelihood  function  has  been  calculated,  the  data  fusion 
process  (association  and  state  estimation)  is  performed  automatically  by  the  Nodestar  engine.  This 
makes  the  Nodestar  engine  extremely  powerful.  It  effectively  reduces  the  problem  of  fusing  the 
information  from  a  new  type  of  sensor  to  a  question  of  modeling  the  physics  of  the  sensor,  a  process 
that  is  well  understood.  Use  of  this  procedure  removes  the  mystery  and  the  need  for  ad  hoc 
procedures  from  data  fusion.  It  provides  a  general  approach  and  framework  for  performing  data 
fusion. 

1.1  Definition  of  Data  Fusion 

To  place  our  discussion  in  a  familiar  context,  we  use  a  modified  version  of  the  definition  of  the 
data  fusion  process  that  has  been  given  by  the  Joint  Directors  of  Laboratories  (JDL)  Data  Fusion 
Subpanel  (see  Waltz  and  Llinas  (1990)).  Data  fusion  consists  of  performing  the  activities  described 
in  the  following  paragraphs. 

Data  association  is  performed  on  the  set  of  reports  (i.e.,  sensor  responses)  produced  by  the 
sensors.  Association  is  the  process  of  partitioning  this  set  into  a  set  of  disjoint  subsets.  Each  subset  in 
this  partition  contains  all  reports  associated  to  a  single  target  and  only  those  reports.  In  some 
situations,  an  additional  subset  is  added  to  contain  the  reports  attributed  to  false  alarms. 

State  estimation  (target  tracking)  is  the  process  of  combining  the  reports  associated  with  a  single 
target  to  produce  an  estimate  of  the  target's  state.  The  target's  state  can  be  very  general,  including 
position,  velocity,  classification,  and  other  target  attributes.  The  estimate  can  be  for  the  target's  past, 
present,  or  future  state.  In  tracking  nomenclature,  these  estimates  correspond  to  smoothing  (past), 
filtering  (present),  and  prediction  (future). 

We  have  combined  attribute  classification,  which  is  listed  as  a  separate  activity  in  the  JDL 
definition,  into  the  state  estimation  process  because,  in  our  view,  attributes  are  properly  part  of  the 
target  state  space. 
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1.2  Related  Activities 

Data  association  and  state  estimation  are  the  core  activities  that  constitute  data  fusion  and  are 
required  for  a  data  fusion  system.  If  one  is  in  a  single  target  situation  (with  no  false  alarms),  then  the 
data  association  part  of  the  problem  may  be  trivial,  but  it  still  takes  place.  The  following  activities  are 
closely  related  to  and  may  be  integrated  into  the  data  fusion  process  in  some  systems.  These  activities 
are,  in  a  sense,  optional  for  a  data  fusion  system. 

Sensor  management  is  the  process  of  using  feedback  from  the  state  estimation  process  (which 
may  include  estimation  of  environmental  parameters)  to  select  sensor  settings  and  modes  of  operation 
that  are  most  effective  for  the  problem. 

Assessment  is  the  process  of  using  the  tactical  picture  produced  by  the  state  estimation  process  to 
provide  estimates  of  the  tactical  situation.  These  estimates  may  include  identification  of  potential 
problems  and  target  intentions. 

Response  planning  uses  the  tactical  picture  and  the  expressed  goals  of  the  decision  maker  to 
recommend  actions.  Allocation  of  sensor  or  processing  resources  to  improve  detection  or 
localization  of  targets  is  an  example  of  response  planning. 

2.  GENERAL  FRAMEWORK  FOR  DATA  FUSION 

In  this  chapter,  we  present  a  general  framework  for  data  association  and  state  estimation.  For  data 
association,  this  framework  will  be  broadly  sketched  with  areas  that  require  further  development 
indicated. 

In  the  case  of  state  estimation,  a  mathematical  theory  based  on  Bayesian  probability  describes 
how  this  should  be  done.  The  principal  impediment  to  the  application  of  this  theory,  heretofore,  has 
been  the  lack  of  computers  with  sufficient  memory  and  speed  to  implement  a  numerical 
approximation  of  the  theory  in  a  satisfactory  manner.  Commonly  available  computers  have  now 
evolved  sufficient  capability  for  the  method  to  be  applied  to  the  full  range  of  tracking  problems. 

2.1  Data  Fusion 

We  have  identified  the  core  activities  of  data  fusion  as  being  data  association  and  state  estimation. 
Before  we  treat  these  as  separate  activities,  let  us  back  up  a  moment  and  view  them  as  part  of  single 
activity,  inference.  Specifically,  the  core  activity  of  data  fusion  is  using  observations  from  various 
sources  to  estimate  the  number  and  state  of  the  targets  in  the  universe  of  interest.  The  observations 
contain  ambiguities,  uncertainties,  and  errors.  The  goal  of  data  fusion  is  to  squeeze  as  much 
information  as  possible  from  these  observations. 

From  the  authors'  perspective,  Bayesian  inference  provides  the  best  framework  in  which  to 
perform  this  inference  task.  The  choice  of  a  Bayesian  framework  is  based  more  on  pragmatic 
grounds  than  on  philosophical  ones.  Of  all  the  common  methods  of  inference  (see,  for  example. 
Pearl  (1988)),  it  provides  the  best  combination  of  a  logical  framework  and  a  powerful  set  of  inference 
tools.  In  addition,  many  of  the  classical  types  of  non-Bayesian  inference  are  equivalent  to  Bayesian 
inference  procedures  with  specific  priors.  This  is  true  for  most  classical  statistical  inference  (Lindley 
1972)  and  for  classical  Kalman  filtering,  too.  Many  of  the  heuristic  inference  systems  have  internal 
inconsistencies  that  make  them  suspect  at  best.  This  is  true  of  most  AI  systems  that  attempt  to  handle 
uncertainties  (see  Pearl  (1988)). 

The  Dempster-Shafer  belief  function  theory  offers  a  logical  approach  to  inference  in  the 
presence  of  uncertainty.  It  is  a  generalization  of  the  Bayesian  probability  model  that  allows 
probability  to  be  assigned  to  a  subset  without  requiring  it  to  be  assigned  to  any  members  of  the 
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subset.  However,  in  its  present  state  of  development,  the  Dempster-Shafer  theory  lacks  the  powerful 
and  complete  set  of  inference  tools  that  Bayesian  theory  provides.  As  an  example,  there  is  not  a  well- 
defined  and  natural  notion  of  expectation  for  belief  functions. 

The  usual  criticisms  of  the  Bayesian  approach  are  that  it  is  computationally  intensive,  and  it 
requires  subjective  judgments  in  terms  of  prior  probability  distributions.  The  cogency  of  the  first 
objection  is  quickly  disappearing  as  the  computational  capability  of  computers  continues  its  rapid 
increase.  The  second  objection  has  substantial  appeal  until  one  considers  the  alternative:  no  prior 
information  is  used  at  all.  One  example,  taken  from  tracking  submarines,  will  illustrate  this  point.  We 
have  often  been  in  situations  where  a  skeptic  will  say,  "I  have  no  prior  knowledge  at  all,  so  how  can  I 
construct  a  prior?"  In  the  case  of  submarine  targets,  there  is  some  obvious  prior  knowledge  that  is 
useful  to  include.  For  example,  the  submarines  are  not  going  to  have  a  speed  faster  than  50  kt,  or 
whatever  the  appropriate  top  speed  is.  They  certainly  are  not  going  to  be  on  land.  Usually,  you 
know  whether  they  are  in  the  Atlantic  or  Pacific  and  in  the  northern  or  southern  hemisphere.  After 
answering  a  few  questions  like  these,  a  rough  prior  can  be  obtained.  Often  a  much  better  prior  can  be 
obtained  with  more  thought.  Using  even  a  rough  prior  is  better  than  having  your  tracker  estimate  a 
submarine's  speed  to  be  100  kt! 

2.2  Bayesian  Inference  Model 

In  very  general  terms,  our  inference  problem  can  be  stated  as  follows.  Let  X  be  the  state  space  of 
our  targets.  The  state  space  can  include  position  (3D),  velocity,  and  any  attribute  information  that  is 
deemed  important,  e.g.,  class  of  ship  or  even  the  ship  name.  We  pick  an  arbitrary  starting  time  and 
designate  it  as  t  =  0.  For  each  t  >,  0,  we  wish  to  know  N(t)  the  number  of  targets  present  and  for  the 
nth  target,  we  want  an  estimate  of  Xfi(t),  the  state  of  that  target  at  time  t.  For  this  discussion,  it  is  more 
convenient  to  think  of  as  a  constant  N,  the  number  of  targets  as  constant  over  the  time  period  of 
interest.  Targets  that  have  not  appeared  in  the  area  of  interest  can  be  assigned  to  a  special  state  until 
they  do. 

Define  the  stochastic  process  5  as  follows: 

S  =  (S(t)  ;  t>0}^{N,X,(t),...,XJt)  ;  t>0}. 

The  prior  information  about  the  targets  and  their  "movements"  through  state  space  is 
incorporated  into  the  probability  measure  defining  this  process.  Thus,  S  becomes  the  prior 
distribution  for  the  Bayesian  inference  process.  Suppose  that  by  time  t,  we  have  a  set  of  observations 
0(t)  =  (Y^,  .  Each  observation  Y^  can  be  a  discrete  event  that  occurred  at  some  time  u 

where  0  <  m  <  t,  or  a  continuous  event  such  as  the  time  series  output  from  a  sensor  that  has  occurred 
over  an  interval  of  time.  The  data  fusion  problem  is  now  equivalent  to  computing  the  posterior 
process 


S  =  {Siu);  u>0\Oit)}  =  {N,X,(u),...,Xj^(u) ;  u>0\O(t)} 

=  {5(m);  u>0}  =  {N,X^iu),...,X^iu);  u>0}. 

The  distribution  of  S(u)=  {N,X^(u),...,X,^(u)  ,u>0}  is  the  posterior  distribution  on  the  state 

of  the  system  at  time  u.  This  is  the  best  estimate  of  S(u)  at  time  u.  The  posterior  distribution 

contains  all  of  the  information  about  5(M)that  we  can  infer  from  0(t).  From  this  distribution  we  can 
compute  means,  86  percent  containment  regions,  approximate  ellipses,  or  whatever  summary 
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information  we  wish  about  the  number  of  targets  and  their  estimated  states.  If  m  =  t,  this  process  is 
called  estimation.  If  u>t,  it  is  called  forecasting,  and  if  m  <  r,  it  is  called  smoothing. 

In  a  sense  we  have  reduced  the  problem  of  data  fusion  to  a  complicated  application  of  Bayes' 
rule.  We  now  know  what  we  would  like  to  compute.  The  question  is:  how  do  we  compute  it?  That  is 
the  topic  of  the  rest  of  this  report. 

2.2.1  Simplifying  Structures 

To  have  any  reasonable  hope  of  computing  the  distribution  of  S ,  we  must  make  some 
simplifying  assumptions  or  approximations.  One  typical  approximation  is  to  split  the  process  of 

computing  S  into  two  steps,  data  association  and  state  estimation.  Let  us  first  consider  the  unified 
activity  for  a  while  longer.  Let  us  assume  that  the  observations  in  0(t)  are  obtained  at  a  set  of  discrete 
times  0<ti<t2<---<tk<  t.  Let 

piSi,...,Si^,s)  =  Pr{S{t^)  =  s^,...,Sit^)  =  s^,Sit)  =  s} 

be  the  prior  probability  (density)  that  the  process  S  passes  through  the  states  si,...,sicat  times  ti,...,tk. 
Suppose  that  the  joint  likelihood  function  L  for  this  set  of  observations  depends  only  on  the  system 
states,  S(ti)^_,,S(tk),  at  these  times.  Then  we  have 

Pr{m)  =  y„...,YM  =  y,\Siu),0<u<t}  =  Pr{Y,(f) 

—  y\j- •  -lYkifk)  —  y^  I *^(^1 ) ~  *^1 ) •  •  ■ ) S{t^. )  —  }. 

We  can  now  define 

L(y„  ...,y,\s„...,s,)  =  Pr{Y,(f)  =  y„ ...,  Y,(t,)  =  y,\S(f)  =  s, . S(t,)  =  sJ_ 

Let  qit,s)  =  Pr{S(t)  =  s}.  The  function  q(t,-)  gives  the  posterior  distribution  on  S(t)  given  0(t). 
By  Bayes  theorem, 

Pr{Y(f)  =  y„...,Y(t,)  =  y,  and S(t)  =  sj 
Pr(Y(f)  =  y„...,Y(tJ  =  yJ 

J ^(y\’  ■  •  ^  ‘^1’  ■  •  •’  ^2’  •  •  •’  ^k’  s}d^\d^2' '  'd^k 

J  ^(yp  •  •  I  •  •)  ^k)p(^\’  ^2’  •  •  •>  ^k’  s)ds^ds2'  "dsi^ds  ^2^ 

Equation  (2)  gives  a  formula  for  calculating  the  joint  probability  (density)  of  S(t)  = 
{N,X^(t),...,XJt)}.  In  general,  there  will  be  no  nice  method  of  computing  Eq.  (2).  Instead,  one 
will  have  to  use  numerical  methods  such  as  Monte  Carlo  integration  to  perform  the  calculations  in 
Eq.  (2).  Observe  that  the  process  of  calculating  q(t,s)  from  Eq.  (2)  does  not  involve  separate  steps  for 
association  and  tracking.  If  the  computations  in  Eq.  (2)  can  be  performed  in  a  feasible  numerical 
manner,  all  of  the  complicated  questions  of  multiple  hypothesis  management  are  subsumed  in  that 
calculation. 
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2.2.2.  Recursive  Computation 

We  now  make  some  additional  simplifying  assumptions  that  will  allow  us  to  perform  the 
calculation  in  Eq.  (2)  in  a  recursive  manner.  First,  the  stochastic  process  S  must  be  Markovian  in  the 
state  space  I+xXx-x  X,  where  /+  is  the  set  of  positive  integers  and  the  product  of  the  X  spaces  is 
taken  a  number  of  times  equal  to  a  reasonable  upper  bound  on  the  number  of  possible  targets. 
Second,  the  likelihood  of  an  observation  Y(ti)  given  the  past  must  depend  only  on  the  present  target 
state,  i.e., 


Pr(Y(tJ  =  y,  I  S(t,)  =  s„ . . .,  S(t,)  =  sj  =  Pr{Y(tJ  =  y,  1  S(tJ  =  sj=  L,(y,  I  s,) , 


and  the  distribution  of  Y{ti)  must  be  independent  of  Y{tj)  given  S(ti)  =  si,...,  S{tk)=  sk  for  i  j. 

Careful  choice  of  the  state  space  is  necessary  to  make  these  assumptions  reasonable. 

Let  po  be  the  probability  (density)  function  for  S(0)  and  pi(si,  si-j)  =  Pr{S(ti)  =  Si  I  S(ti.i)  =  j/. 
;}.  Then 


p(s^,...,Si^)  =  >Si.i)Po(So)dsQ, 

i  =  \ 

k 

L(y„ . . ., yj 5i, .. ., 5 J  =  YihiYi ' 


Substituting  Eq.  (3)  into  Eq.  (2),  we  obtain 


(3) 


2  k-\  k-\ 

(l(tk>s}  =  —Ah(yk\  ^  A  ( Yi  I  )Pk  i  s,  )]!  Pi  ( ^i>  ^i-i  ^Po  ( ^0  )dSods,  ■  •  •  ds,_2ds,_, 

c  i=i  i=j 


1  r  f 

=  —  4  f y,  I  ^ jj  p,(s,s,_J  J  n  A ( Yi  1  )Pi  ( )Po ( ^0 )dSodSy  ■  •  ■  ds,_2 

=  —  L^fy*  I s)j Pj^(s,Si^_^)q(ti^_^, , 


d^k-\ 


(4) 


where 


C  =  \\h(Yk\  Pk(s>  h-M  h-v  Sk-x)dSk-x 


ds. 


Equation  (4)  provides  a  recursive  method  of  computing  q(tk,s).  The  integral  on  the  right-hand 
side  of  the  last  equality  in  Eq.  (4)  performs  the  "motion  update"  from  the  time  tk-i  of  the  last 
computation  to  the  present  time  tk-  The  multiplication  by  Lk(yk^s)  performs  the  information  update. 
Division  by  C  renormalizes  the  numerator  to  a  probability  distribution.  If  there  has  been  no 
observation  at  time  tk,  then  Lk(yk^s)  =  1,  and  there  is  no  information  update,  only  a  motion  update. 

2.3  Separation  of  Data  Association  and  State  Estimation 

Most  data  fusion  systems  separate  the  activities  of  data  association  and  state  estimation.  From  this 
point  on,  this  report  returns  to  that  perspective.  The  general  approach  is  to  first  perform  the 
association  by  developing  some  method  of  partitioning  the  set  of  observations  into  a  set  of  disjoint 
subsets  that  cover  the  entire  set  of  observations.  The  observations  in  a  subset  are  associated  to  a  single 
target.  No  other  observations  are  associated  with  that  target.  Having  performed  the  associations,  the 
data  fusion  problem  is  reduced  to  a  set  of  estimation  or  tracking  problems,  one  for  each  subset. 
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The  process  usually  proceeds  recursively.  It  begins  with  a  set  of  observations.  These  are 
associated  to  form  targets  and  state  estimates  (tracks)  are  produced.  When  the  next  set  of 
observations  is  obtained,  these  tracks  are  updated  to  the  time  of  the  observations  to  provide 
information  to  the  association  process.  The  choice  of  how  to  associate  the  observations  to  existing  or 
new  tracks  will  often  be  ambiguous.  Most  systems  require  that  a  single  choice  be  made.  Others 
consider  multiple  choices  (partitions  of  the  set  of  observations)  and  assign  a  probability  to  each 
partition  or  hypothesis.  These  systems  are  called  multiple-hypotheses  correlator-trackers.  The  term 
multiple- hypotheses  refers  to  association  hypotheses.  Even  though  the  processes  of  association  and 
state  estimation  take  place  in  separate  steps,  they  are  very  closely  related.  In  particular,  the  motion- 
updated  tracker  estimates  will  affect  the  association  of  observations  to  tracks,  and  obviously  the 
association  of  observations  to  tracks  affects  the  state  estimates. 

The  paper  by  Mori  et  al.  (1986)  describes  a  Bayesian  framework  for  a  two-stage  process  of 
association  and  estimation  that  provides  for  multiple  association  hypotheses.  Given  that  one  has 
chosen  to  view  the  association  and  estimation  steps  as  separate  although  related  activities,  this  paper 
provides  an  excellent  approach  to  data  association  that  is  far  more  powerful  and  flexible  than  that 
provided  by  the  Kalman  filter-based  methods.  The  interested  reader  is  referred  to  this  paper  for  the 
details  of  the  approach.  Although  Mori  et  al.  do  discuss  nonlinear  state  estimation  procedures,  the 
discussion  and  examples  given  in  the  paper  do  not  show  the  full  power  of  this  approach  nor  do  they 
provide  a  general  methodology  for  implementing  a  nonlinear  state  estimation  (tracking)  process. 

The  remainder  of  this  report  concentrates  on  describing  an  approach  to  nonlinear  state  estimation 
and  describing  a  successful  and  powerful  engine.  Nodestar,  that  implements  that  approach. 

3.  BAYESIAN  FILTERING 

We  now  turn  to  the  problem  of  state  estimation.  Following  the  separation  principle  described  in 
the  previous  section,  we  take  the  view  that  the  observations  have  been  associated,  i.e,.  partitioned  into 
subsets  with  each  subset  corresponding  to  a  single  target.  The  state  estimation  problem  is  to  use  the 
observations  in  a  given  set  to  estimate  the  state  of  the  target.  We  wish  to  make  a  distinction  at  this 
point.  We  call  targets  the  actual  objects  whose  state  we  are  estimating.  For  example,  in  a  surface  ship 
tracking  problem,  the  targets  are  the  ships.  A  track,  on  the  other  hand,  is  the  set  of  observations  that 
are  associated  into  a  single  subset  under  the  partitioning  performed  by  the  association  algorithm.  The 
distinction  is  that  targets  are  real  objects.  Tracks  are  constructs  created  by  the  data  fusion  system. 

In  this  section,  we  describe  a  state  estimation  procedure  called  Bayesian  filtering.  In  a  sense, 
Bayesian  filtering  is  not  a  new  concept.  It  is  very  closely  related  to  probabilistic  filtering  as  defined 
by  Jazwinski  (1970).  This  procedure  is  often  mentioned  in  standard  texts  and  papers  on  data  fusion, 
such  as  Bar-Shalom  and  Fortman  (1988),  Mori  et  al.  (1986),  and  Waltz  and  Llinas  (1990).  However, 
it  is  quickly  dismissed  as  impossible  to  implement,  and  the  references  move  on  to  extended  Kalman 
filtering  or  weighted  Gaussian  sum  approaches.  One  purpose  of  this  discussion  is  to  reintroduce  this 
approach  and  to  say  that  it  is  feasible  and  has  been  implemented  in  a  working  system.  Nodestar. 
Nodestar  has  been  tested  in  operational  situations  on  surveillance  problems.  Furthermore,  Nodestar  is 
not  simply  a  system  designed  for  a  single  problem.  It  is  an  engine  that  can  readily  be  modified  to 
handle  a  wide  range  of  tracking  problems  with  diverse  data  sources. 

The  fact  that  Nodestar  demonstrates  the  computational  feasibility  of  implementing  a  Bayesian 
filter  radically  alters  the  nature  of  the  data  fusion  problem  faced  by  the  R&D  community.  In  the 
past,  the  community  has  spent  large  amounts  of  effort  trying  to  adapt  Kalman  filtering  to  handle 
disparate,  nonlinear  types  of  information.  The  inherent  limitations  of  Kalman  filtering  doom  this 
effort  to  failure,  but  researchers  have  continued  to  pursue  this  approach  because  they  have  felt  that 
Kalman  filtering  provided  the  only  computationally  feasible  approach.  Nodestar  has  changed  all 
this.  Bayesian  filters  have  been  shown  to  be  computationally  feasible  and  extremely  flexible  in 
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fusing  diverse  types  of  information.  Furthermore,  the  fusion  process  proceeds  according  to  a 
consistent  and  correct  mathematical  model.  In  view  of  this,  we  feel  that  researchers  should  redirect 
their  efforts  toward  using  the  Bayesian  filtering  methodology  to  develop  systems  that  handle  these 
data  fusion  problems  in  a  natural,  straightforward,  and  mathematically  consistent  fashion. 

The  use  of  Kalman  filtering  should  be  limited  to  high  data  rate  situations  in  which  the  linear- 
Gaussian  assumptions  of  Kalman  filtering  are  reasonably  satisfied.  Tracking  aircraft  with  radar  is  an 
example  of  such  a  situation. 

3.1  Formulation 

Bayesian  filtering  is  based  on  the  mathematical  theory  of  probabilistic  filtering  described  by 
Jazwinski  (1970).  Bayesian  filtering  is  a  method  of  estimating  a  target's  state  at  a  sequence  of  times. 
The  target's  state  vector  typically  includes  position  but  may  also  include  velocity  and  other  variables 
such  as  those  representing  target  attributes  or  characteristics.  Let  X  be  the  target's  state  space.  The  a 
priori  information  about  the  target  is  represented  by  a  stochastic  process  {X(r);  t  >  0},  where  X{t)  is 
the  target  state  at  time  t.  Sample  paths  of  this  process  correspond  to  possible  target  paths  through  the 
state  space  X.  For  any  discrete  set  of  times  0  <  ti  <  r2  -  •  •<  t/j-i  <  we  define 

Xn)  =  Pr{Xiti)  =  XI,...,  X{tn)  =  (5) 

to  be  the  prior  probability  (density)  of  the  sample  paths  of  the  process  that  pass  through  the  states 
jci,...,  Xfi  at  times  ti,...,tn-  Suppose  that  we  obtain  observations  0(r„)  =  {y(ti),...,  F(r„)}  at  times 
ti,...,tfi.  In  general,  these  observations  are  noisy,  nonlinear  functions  of  the  target  state  or  path.  We 
assume  that  we  can  calculate  the  joint  likelihood  L  for  these  observations  conditioned  on  target  path 
and  that  L  depends  only  on  the  target  state  at  times  t ],..., tn-  Thus, 

Pr{0(tn)=  (yi,-..,  Jn)  I  X(u),  0<u<tn}  =  Pr{0(tn)  =  (yi,---,  yn)  I  X(ti),...,  X(tn)}. 


Let 


L(yi,...,  yn\xi,...,  %)  =  Pr{G(t„)  =  (yi,...,  y„)  I  Z(fi)  =  xi,...,  X(t„)  =  x„}. 


(6) 


The  objective  of  Bayesian  filtering  is  to  calculate 

q{tn,  x)  =  Pr{X{tn)  =  x  I  0{tn)  “  (yi,...,  yn)}  for  X  e  X,  (7) 

the  posterior  distribution  on  target  state  at  time  t^  given  the  information  0(fn).  The  function  q{tn,  ■)  is 
the  complete  solution  to  the  filtering  problem  in  the  sense  that  it  embodies  all  the  information  about 
X{tn)  that  is  contained  in  the  observations  0{tn)  and  the  prior  distribution  on  {X(0;  t  >  0}.  Using 
Bayes'  theorem,  we  obtain 


Pr(0(tJ  =  {y„...,yJmdX(tJ^x} 
Pr{0(tJ^(y„....yJ} 


f  L(  y,, . . .,  yj  Xp . . x„_i,  x)p(  x,,  x^ . . .,  x„_i,  x)dx^dx.^--  •c(x„_i 

=  - - 5 - for  xeS. 

J  L( yp . . ., y„  I  Xp . . .,  xjp( Xp  Xj, . . ., xjdx^dx^- • 


(8) 
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In  most  cases,  we  are  forced  to  compute  q(tn,  x)  numerically.  The  following  subsections  describe 
the  most  common  ways  of  doing  this. 

3.1.1  Batch  Computation 

In  the  batch  method,  one  computes  the  integrals  in  Eq.  (8)  each  time  a  new  observation  is 
obtained.  Most  often  a  Monte  Carlo  evaluation  of  these  integrals  is  performed  by  replacing  the 
stochastic  process  {X(r);  f  >  0}  by  a  finite  number  M  of  sample  paths  drawn  from  the  process. 
Typically  M  is  about  1000.  Usually  the  paths  x„)  are  chosen  so  that  x„)  =  1/M. 

Having  done  this,  one  computes  the  likelihood  function  L  for  each  of  these  paths  as  a  function  of  the 
observations,  0(r„).  The  numerator  in  Eq.  (8)  is  the  sum  of  these  likelihood  functions  over  all  the 
paths  that  end  at  x  multiplied  by  the  constant  factor  1/M.  The  denominator  is  the  same  sum  taken 
over  all  paths.  Batch  computation  is  the  only  method  of  computing  Eq.  (8)  at  this  level  of  generality. 

3.1.2  Recursive  Computation 

When  possible,  it  is  more  efficient  to  compute  q{tn,  x)  in  a  recursive  rather  than  batch  fashion.  To 
do  this,  we  require  the  following  additional  assumptions.  First,  the  stochastic  process  {X(0;  t  ^  0} 
must  be  Markovian  in  the  state  space  X.  Second,  the  likelihood  of  the  observation  Y{ti)  given  the  past 
must  depend  only  on  the  present  target  state,  i.e., 

Pr{  Y(ti)  =  yi  I  Z(ti)  =  XI ,. . . ,  =  x/}  =  Pr{  Y{ti)  =  yt  I  X( tO  =  x;  }  (9) 


=  Liiyi  I  Xi), 


and  the  distribution  of  Y{ti)  must  be  independent  of  Y{tj)  given  X{t\)  -xi,...,  X{tn)-  Xn  for  i^j. 
Careful  choice  of  the  state  space  usually  makes  these  assumptions  reasonable. 

Let  pQ  be  the  probability  (density)  function  for  X(0)  and  pfxi,  xpi)  =  Pr{X{ti)  =  x/ 1  X{ti.\)  = 


xpi).  Then, 

p(x^,...,xj  =  Xi-i)Po(Xo)dXo, 

n 

A  { yi  I  )• 

i~\ 


(10) 


Substituting  Eq.  (10)  into  Eq.  (8),  we  have 


1  Ai-i 

^(^n>  X)  =  —  J  4  (')'«  I  A  (  yi  '  ( X,  X„_i  jfl  Pi  ( ^i’  ^i~.\ )  Po  ( ^0  jdx^dx^  •  *  •  dx^_2dx^ 

C  ^  i=\  i=i 

1  r 

=  — I  x)j Pjx,  x„_,  j  jn[44,'  I  x,)Pi(x„  X._,)]p,( x,)dx,dxp- 

^  i=\ 


dx, 


n-l 


^4(y„|x)J p„{x,x„_,)q(t„_„x„_,)dx„_,  for  x  e  X., 


(11) 


C  =  Pn(x,x„_i)q{t„_„x,^_^)dx,,_^ 


dx. 


where 
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The  recursion  in  Eq.  (11)  breaks  into  two  steps,  motion  updating  and  information  updating. 

3.1.3  Motion  Updating 

The  integral  in  the  last  line  of  Eq.  (11)  updates  the  target  state  distribution  at  time  tn-\  for  the 
"motion"  of  the  target  through  the  state  space  up  to  time  Let  q'{tn,-)  denote  the  motion-updated 
distribution.  Then, 

=  \  forXG^  (12) 

3.1.4  Information  Updating 

The  pointwise  multiplication  of  q'(tn,-)  by  the  likelihood  function  L„(y„,  ■)  and  the  normalization 
by  C  constitutes  the  information  update.  Thus 

=  ^k(ynV)(f(K^x)  forx  e  X  (13) 


Observe  that  sensor  information  such  as  detections,  no  detections,  elliptical  position  estimates, 
estimates  of  the  bearing  of  the  target  from  a  sensor,  measured  signal-to-noise  ratio,  observed 
frequency  of  signal,  and  many  more  can  be  integrated  into  the  tracking  solution  via  the  likelihood 
functions  L^.  If  there  are  several  types  of  information  or  measurements  received  at  a  single  time,  then 
Ln  becomes  a  joint  likelihood  function  for  all  the  information  received. 

3.1.5  Nonlinear  Tracking 

Equations  (12)  and  (13)  do  not  require  the  observations  to  be  linear  functions  of  the  target  state. 
They  do  not  require  the  measurement  errors  or  the  probability  distributions  on  target  state  to  be 
Gaussian.  We  call  this  nonlinear  filtering  or  tracking. 

3.1.6  Monte  Carlo  and  Markov  Chain  Implementations 

At  this  level  of  generality,  Eqs.  (12)  and  (13)  must  be  evaluated  numerically.  The  existence  of 
high-powered  scientific  workstations  allows  us  to  compute  and  display  tracking  solutions  for  complex 
nonlinear  trackers.  There  are  two  standard  approaches  to  numerical  computation  of  these  equations. 
In  the  Monte  Carlo  approach,  the  stochastic  process  {X(t);  r  >  0}  is  approximated  by  recursively 
simulating  a  number  M  of  target  paths  through  the  state  space  S.  Typically,  each  path  is  assigned  an 
a  priori  probability  1/M.  At  each  time  step  t^,  this  path  is  moved  forward  to  its  position  Xn  at  time 
and  its  probability  is  multiplied  by  L(y«  I  Xn)IC.  Moving  the  paths  forward  to  tn  produces  a  discrete 
numerical  approximation  to  Eq.  (12).  Multiplication  by  L{yn  I  Xn)/C  accomplishes  Eq.  (13). 

In  a  Markov  chain  approach,  one  typically  discretizes  the  state  space  so  that  Eq.  (12)  is  computed 
through  the  use  of  discrete  transition  probabilities.  The  likelihood  functions  are  also  computed  on 
the  discrete  state  space.  Markov  chain  trackers  often  operate  with  more  than  300,000  states.  To 
cover  a  wide  range  of  possible  situations,  they  use  a  dynamic  grid  for  the  state  space.  The  size  of  the 
cells  in  the  grid  adjusts  as  the  distribution  on  target  state  expands  or  contracts. 

3.1.7  Kalman  Filters 

A  Kalman  filter  tracker  represents  a  special  case  of  the  recursion  in  Eqs.  (12)  and  (13)  that  holds 
under  special  assumptions.  The  target  motion  process  {X(r);  f  ^  0}  must  satisfy  a  linear  stochastic 
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differential  (or  difference)  equation,  and  the  distribution  of  X(0)  must  be  Gaussian.  The  observations 
Y{ti)  must  be  linear  with  Gaussian  errors,  i.e., 

Y{ti)  =  M{ti)xti  +  e,  for  /  =  1 , . . . ,  n,  (14) 

where  Mitt)  is  a  matrix  and  the  e/  are  independent,  mean  0,  Gaussian  random  variables  for  i  = 

Under  these  assumptions,  Jazwinski  (1970,  Chapter  7)  shows  that  the  recursion  in  Eqs.  (12)  and  (13) 
reduces  to  the  standard  recursive  Kalman  filter.  This  is  called  a  linear-Gaussian  filter. 

The  linear-Gaussian  restrictions  are  severe  in  environments  where  one  is  required  to  fuse  diverse 
types  of  information.  These  often  involve  non-Gaussian  distributions  and  nonlinear  relationships 
between  the  measurements  and  the  target  state.  A  line  of  bearing  measurement  is  an  example;  so  is  a 
detection/no  detection  observation.  Many  Kalman  filters  handle  this  problem  by  approximating 
nonlinear  measurement  relations  with  linear  ones.  This  is  called  an  extended  Kalman  filter.  The 
success  of  this  approximation  varies.  Linear  motion  models  are  also  restrictive.  For  example,  one 
cannot  have  a  velocity  distribution  that  varies  over  space. 

3.2  Bayesian  Prediction  and  Smoothing 

In  filtering  terminology  (Jazwinski  {\91Q)),  filtering  is  the  process  of  estimating  the  present  target 
state,  prediction  is  estimating  the  target's  state  at  some  time  in  the  future,  and  smoothing  is  estimating 
the  target's  state  for  a  time  in  the  past.  In  this  section,  we  discuss  the  extension  of  the  above  methods 
to  the  problems  of  prediction  and  smoothing. 

In  a  sense,  the  prediction  problem  is  a  trivial  extension  of  the  filtering  solution  that  we  have 
already  obtained.  One  simply  applies  the  motion  update  step  given  in  Eq.  (12)  to  the  funcdon  q(t,-), 
which  gives  the  posterior  distribution  on  target  state  at  the  present  time  t,  to  compute  the  distribution 
on  target  state  at  any  time  u>t.  For  convenience,  we  consider  situations  in  which  there  is  information 
about  the  target's  future  location  that  we  have  not  included  in  the  present  state  estimate.  In  this  case, 
the  problem  of  prediction  becomes  more  interesting  and  is  parallel  to  the  problem  of  smoothing. 

3.2.1  Formulation 

We  limit  ourselves  in  this  discussion  to  the  situation  where  the  computation  can  be  performed 
recursively.  In  practice,  this  is  the  only  situation  in  which  it  is  feasible  to  perform  the  calculations 
involved. 

Let  [0,T]  be  the  time  interval  of  interest,  and  let  0  <  tj  <  t2...<  tn-j  <  tn^The  a  discrete  set  of 
times  within  this  interval.  Usually,  these  are  times  at  which  observations  have  taken  place.  We  assume 
the  process  {X(r);  r  >  0}  is  Markovian  and  that  0(T)  ={y(ri),...,  Y(tn)}  is  the  set  of  observations 
obtained  in  [0,T].  Paralleling  Eq.  (9),  we  assume  that 

PrfY(ti)  =  yi  I  X(ti)  =  XI,...,  X(tn)  =  xj  =  Pr{Y(ti)  =  yt  I  X(ti)  =  xd  =  Li(yi  I  xO  (15) 

and  that  Eq.  (10)  holds. 

For  1  <k  <n,  define 


q[t„xp[T)) 


Pr\ 

[o(r) : 

=  (yi.- 

..,y„)andX(4)  =  x} 

Pr 

{0(7-) 

} 

(16) 


Technical  Documentation  of  Nodestar 


11 


Note  that  q(t„, x)  =  q(t„,x\0( t„ )) .  For  some  constant  C" , 

q(t^ ,  |0(r))  =  ■—  J  [^(^I  \^i  )Pi  >  ^i-I  )]Po  (^0  )^o  ' '  ■  ■ 

V  1=1 

Observe  that  we  do  not  integrate  over  xk  on  the  right-hand  side  of  Eq.  (17).  We  can  write 
q(t^,x,^\0(T)) 

=  -^L(y^\xjj  Pk(x^,x,jYl[L(y,\Xi)p,(x,,x,_,)]pJXo)dXf,---dx^_, 

C  i=i 

^ltl[Uy,  \Xl)Pi(Xi,Xl_^)]dx^^^■^■dx^ 

i=k+\ 

=  -^L(y,\x^  )q-(t„  X, )J  )',■ '  h)Pi( ^i>  h-i 

^  i~k+\ 

where  q‘  is  defined  in  Eq.  (12)  and  C"  is  the  appropriate  normalizing  factor. 


Mxk\yk^v-yn)  =  J  tl[L{yi\xi)Pi(xix,_,)]dx,,,-dx„ . 

(=Jt+l 

Then 


(17) 


fk(^k  ~  ^k  ^  ^(h+l)  ~  Tjl:+1<'”>  ^(^n)  ynl  ^ 

We  can  now  write 

Observe  that  q  is  proportional  to  a  product  of  factors  corresponding  to  the  effect  of  past 
information,  q~(t^,xj\  the  future  information,  I (^t+p •••>'„)) ;  and  the  present  information 

L(y^\xJ.  Specifically,  q~(t^,xj  is  the  probability  density  that  X(tic)  =  xk  conditioned  on  the  past 
information  O(tk-i).  Note  that  q~(t^,-)  includes  the  contact  information  for  times  t  <  tk  and  is 
"motion-updated"  from  tk-i  to  tk-  Uyb')  ‘s  the  likelihood  function  representing  the  information  in 
the  observation  Y{tk)  =  yk  at  the  present  time,  and  y^^,,...yj  is  proportional  to  the 

probability  density  that  X(tk)  =  xk  given  the  future  observations  .■.,Y(tJ  =  y„.  We  call 

fk  the  future  likelihood  function. 

We  have  already  shown  how  to  compute  q''(t^,x^)  recursively.  There  is  a  backward  recursion 
that  one  can  use  to  compute  . .  .y„) .  In  particular,  for  x  e  X , 
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/„(:^|0)  =  1 

fn-i(x\yn)  =  jL(yJxJp^ix^,x)dx^  (19) 

A  I  >'*+1  -  -  )  =  j  A+1  (x  lyi+2  V  • ,  I  )Pm  (Xm  ’  x)dx,^i 

fotO<k<  n-2. 


3.2.2  Smoothing 


The  function  qit,-\0{tj)  for  0<t<T  gives  the  probability  distribution  representing  the 
smoothed  estimate  of  the  target's  position  at  time  t  given  the  observations  0(T).  Let 

ll(f)=  [  xq(t,x\0(T))dx  for  0<t<T .  Then  |X(f)  is  the  mean  of  the  smoothed  estimate  of  the 
Jx 

target's  position  at  time  t.  One  could  plot  p  over  [0,7]  to  present  an  estimated  smoothed  path  for  the 
target  given  the  complete  set  of  observations  in  0(7). 

In  the  special  case  where  the  Kalman  filtering  assumptions  hold,  the  function  q~(t^,-)  is  the 
normal  probability  density  that  results  from  running  a  forward  Kalman  filter  on  the  observations 
O(tic-i)  and  motion-updating  to  time  t^.  L(y^l-)  is  a  normal  density  function  corresponding  to  the 
observation  at  tk,  and  proportional  to  the  output  of  a  backward  Kalman  filter 

operating  on  the  contacts  K(4^,)  =  F(f„)  =  One  important  observation  about  the 

presentation  given  above  is  that  it  shows  how  to  compute  smoothed  tracks  for  nonlinear  filters. 
Smoothed  tracks  are  usually  of  interest  when  one  is  reconstructing  target  tracks  after  the  fact. 

3.2.3  Land  Avoidance 

One  version  of  the  Nodestar  system  described  below  is  designed  for  tracking  submarines.  We 
know  for  certain  that  submarines  will  not  travel  on  land.  They  go  around  islands  and  avoid  coast 
lines.  Most  tracking  programs  do  not  account  for  this  piece  of  information.  Kalman  filters,  for 
example,  are  quite  capable  of  projecting  submarine  tracks  across  land  and  even  continents.  In  the 
submarine  surveillance  version  of  Nodestar,  we  account  for  this  behavior  by  developing  a  set  of  land 
avoidance  likelihood  functions  that  modify  the  probabilities  on  the  position  and  velocity  components 
of  the  target's  state  as  it  approaches  land.  The  effect  of  this  is  to  modify  the  motion  model  used  by 
the  system.  The  advantage  to  doing  this  by  likelihood  functions  is  that  these  functions  can  be  applied 
to  any  underlying  motion  model.  This  means  that  we  do  not  have  to  reconstruct  the  motion  models 
to  make  sure  they  avoid  land  each  time  we  change  the  underlying  model.  It  will  happen 
automatically  by  use  of  the  land  avoidance  likelihoods. 

The  formalism  is  entirely  parallel  to  that  which  leads  to  Eq.  (18).  For  each  time  t^,  the  likelihood 
function  L( I  •)  is  the  land  avoidance  likelihood  function.  Since  this  function  is  the  same  for  all 
times,  let  us  call  this  function  La-  The  value  of  this  function  is  small  for  those  states  close  to  land  and 
for  those  velocities  that  are  likely  to  put  the  target  close  to  land  before  a  velocity  change  takes  place. 
With  this  in  mind,  Eq.  (18)  provides  a  way  of  computing  the  distribution  on  target  location  at  any 
time  tk  knowing  that  it  avoided  land  in  the  past,  is  avoiding  it  in  the  present,  and  will  avoid  it  in  the 
future.  Since  the  land  avoidance  is  known  ahead  of  time,  we  are  actually  performing  a  prediction 
even  though  the  formalism  is  that  of  smoothing. 

If  an  observation  is  obtained  from  a  sensor  at  any  time  tk,  its  likelihood  function  is  calculated  and 
combined  by  pointwise  multiplication  with  the  land  likelihood  function  to  form  a  joint  likelihood 
function  Lj^  to  be  applied  at  time  t^.  With  this  addition,  the  computations  leading  to  Eq.  (18)  tire 
unchanged. 
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Often,  the  computational  load  involved  in  calculating  the  factor  I •••)'„)  incorporating 

the  land  likelihoods  in  the  future  is  too  burdensome.  In  this  case,  we  use  the  recursion  in  Eqs.  (12) 
and  (13)  but  with  each  observation  likelihood  Ln  multiplied  pointwise  by  the  land  avoidance 

likelihood  function. 


3.2.4  Modified  Motion  Model 

The  application  of  a  land  likelihood  function  La  to  condition  the  future  and  present  motion  of 
the  target  is  equivalent  to  modifying  the  transition  probabilities  of  the  Markov  motion  model.  To  see 
this,  let  the  future  likelihood  function /„  be  defined  by  Eq.  (19)  for  /:  =  0,...,  n  with  l-j  =  La  for 
k  =  0,...,  n.  Define  a  modified  initial  probability  measure  pi  and  modified  transition  probabilities 
p*  for  /  =  1,...,  n  as  follows: 


pI(Xo)  = 


Po(^o)fo(^o)^a(^o) 

j^Po(x')fo(x')LJx')dx' 


for  Xq  e  X, 


=  fotx,_„x,eX. 


(20) 


Then,  for  any  path  (xq . x„),  we  have 

pl(Xo)nPi^’  x._,)LJx,) 

Jof^oy  1=1 


where 


1  ^ 

“  Po( (^0  Pi (^i* 

^  1=1 


C  —  J  Po(  ^o)f o(  ^o)^a( 

X 

=  Ipo(xJLJxJ  \ 

X  L  «=i 

Xi.Ofi(Xi)LJx^)dXo---dx„. 


\dXn 


(21) 


X  X  <=0 


(22) 


From  Eqs.  (21)  and  (22),  we  can  see  that  the  probability  measure  on  the  path  (xq,...,  x„)  produced  by 
the  modified  transition  probabilities,  p*,  is  exactly  the  conditional  probability  of  the  path  (xq,...,  x^) 
given  the  land  likelihood  at  times  tg,...,  tn- 

4.  NODESTAR  ENGINE 

Nodestar  is  a  discrete,  nonlinear,  non-Gaussian,  Bayesian  filter.  It  can  process  disparate  forms  of 
information  on  multiple  targets  at  speeds  sufficient  for  surface  ship  and  submarine  applications. 
Nodestar  implements  the  theory  by  representing  distributions  and  likelihood  functions  as  values  on  a 
discrete  grid  with  dynamically  changing  resolution.  Currently  available  computers  represent  most 
forms  of  information  applicable  to  surface  ship  and  submarine  tracking  with  appropriate  accuracy. 
The  computations  required  by  the  theory  can  be  performed  with  sufficient  speed  to  allow  real-time 
processing. 
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In  the  following  sections,  we  describe  the  Nodestar  engine  in  general  and  the  Spotlight 
implementation  in  particular.  Even  though  the  emphasis  in  the  following  discussion  is  on  the 
Spotlight  implementation,  the  reader  should  note  that  the  engine  itself  is  general  purpose  in  that  the 
motion  models  and  likelihood  functions  described  in  this  section  and  in  Section  5  can  be  easily 
replaced  by  ones  that  are  appropriate  for  markedly  different  applications.  We  provide  a  hint  of  that 
capability  in  Section  7. 

4.1  State  Space 

Nodestar  uses  a  discrete  state  space  of  cells  with  dynamic  resolution.  The  typical  state  space  has 
six  or  more  dimensions  and  300,000  to  500,000  cells.  Each  cell  contains  a  probability.  The 
collection  of  probabilities  in  the  cells  is  the  joint  probability  distribution  for  the  target  state  at  a  given 
time.  In  the  Spotlight  version  of  Nodestar,  the  state  space  is  six-dimensional: 

Position  (3D):  The  dimensions  are  latitude,  longitude,  and  depth.  The  system  operates  in 
spherical  coordinates  so  that  the  latitude-longitude  cells  at  a  given  time  have  a  fixed  width,  Aion 
and  height  Aiat  in  latitude  and  longitude.  This  means  that  the  cells  will  be  somewhat  narrower 
in  the  northern  parts  of  the  grid  than  in  the  southern  parts.  The  motion  models  and  likelihood 
functions  in  Nodestar  are  based  on  spherical  geometry  so  that  all  calculations  are  performed 
on  the  sphere.  Displays  can  be  presented  in  any  projection  desired.  The  Spotlight  version  of 
Nodestar  uses  a  Mercator  projection. 

Velocity  (2D):  The  spotlight  version  of  Nodestar  uses  a  2D  velocity.  This  is  satisfactory  for 
submarine  or  surface  ship  problems.  The  velocity  is  the  local  velocity  along  a  rhumb  line 
course.  Appendix  A  defines  a  rhumb  line. 

Attribute:  Spotlight  uses  the  class  of  the  submarine  target  as  an  attribute.  This  attribute  along 
with  a  database  specifying  the  acoustic  and  radar  characteristics  of  each  class  allows  Nodestar 
to  calculate  many  of  the  likelihood  functions  described  below. 

One  could  substitute  the  ship  name  or  any  other  attribute  for  target  class  if  the  appropriate 
database  is  available  to  support  it.  Other  attributes  can  be  added.  Each  additional  attribute  increases 
the  dimension  of  the  state  space,  but  this  is  a  computational  rather  than  conceptual  problem. 

4.2  Dynamic  Resolution 

The  resolution  of  the  lat-lon  and  velocity  grids  changes  dynamically  as  the  size  of  the  probability 
distribution  represented  by  the  grid  changes.  These  dynamic  grids  are  represented  by  the  Nodestar 
window  and  velocity  wheel. 

4.2.1  Spatial  Regridding 

At  each  time  t.  Nodestar  calculates  the  marginal  two-dimensional  distribution  on  the  lat-lon  cells. 
For  each  lat-lon  cell,  it  sums  the  probability  from  all  states  whose  lat-lon  value  falls  in  that  cell.  For 
each  track  (target),  there  are  a  fixed  number  of  lat-lon  cells,  typically  50  x  50.  The  location  of  these 
cells  is  determined  by  a  center  (lat.  Ion)  and  the  values  of  Aion  and  Aiat-  The  probabilities  in  these 
cells  are  displayed  as  color-coded  on  a  monitor.  The  display  takes  place  in  the  Nodestar  window  for 
that  track.  This  window  is  superimposed  on  a  geographical  display.  Usually,  Ajon  4^  Aiat  so  that  the 
window  is  rectangular  rather  than  square.  When  regridding,  Nodestar  can  change  the  length  of  the 
cell  independently  of  the  width.  This  is  useful  if  the  probability  distribution  is  not  square.  For 
example,  consider  the  probability  map  in  Fig.  1. 
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Fig.  1—  Grid  for  spatial  probability  map 

Suppose  that  we  would  like  to  expand  the  coverage  of  the  probability  map  in  the  east-west  but 
not  in  the  north-south  direction.  For  example,  suppose  that  we  would  like  to  change  our  cells  so  that 
they  look  like  Fig.  2.  In  this  example,  we  redefined  the  cells  by  doubling  the  width  and  keeping  the 
height  constant.  We  call  this  regridding  by  the  factors  (2.0,  1.0).  With  Nodestar,  any  pair  of  positive 
numbers  may  be  used  for  regridding  factors. 


Pig^  2 —  Grid  of  new  cells 


The  size  and  the  center  of  the  window  are  determined  dynamically  by  the  location  and  spread  of 
the  marginal  distribution  on  lat-lon  in  the  manner  described  below.  The  result  is  that  the  center  of 
this  window  shifts  as  the  center  of  distribution  moves.  This  size  of  the  window  will  shrink  as  the 
probability  distribution  becomes  more  concentrated,  but  the  number  of  lat-lon  cells  remains  constant. 
Thus,  as  the  target  distribution  becomes  more  localized,  the  size  of  the  grid  cells  becomes  small  and 
the  resolution  of  the  grid  becomes  finer.  As  the  distribution  spreads  out,  the  reverse  happens. 

The  notion  of  splitting  motion  updates  into  steps  was  developed  in  response  to  a  problem  that 
Nodestar  encounters  when  the  probability  distribution  is  small  or  the  motion  update  time  interval  is 
long.  Unless  care  is  taken  in  these  cases,  the  translation  portion  of  the  motion  update  moves  a 
significant  part  of  the  probability  mass  (or,  in  extreme  cases,  all  of  the  mass)  outside  the  Nodestar 
grid,  and  this  probability  is  lost. 

The  solution  to  this  problem  is  to  break  the  update  interval  into  steps  that  allow  Nodestar  to 
expand  the  size  of  the  distribution  by  a  specific  fraction  per  step.  At  the  end  of  each  step.  Nodestar 
regrids  the  distribution  so  as  to  create  a  buffer  of  empty  cells  around  the  distribution  that  is  large 
enough  to  accommodate  the  next  motion  step.  The  accuracy  and  speed  of  the  new  method  depends 
heavily  on  the  fraction  by  which  the  distribution  is  allowed  to  grow.  This  fraction  is  one  of 
Nodestar’s  configuration  parameters,  the  Single  Step  Motion  Threshold.  Legal  values  are  between 
0.01  and  0.4. 
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The  advantage  to  the  above  algorithm  is  that  no  mass  will  be  lost  because  of  motion  updates.  Of 
course,  one  does  not  necessarily  want  to  keep  the  lowest  probability  events  inside  the  Nodestaf  grid, 
so  the  spatial  regridding  algorithm  also  has  a  parameter  that  controls  how  much  mass  it  may  discard 
before  creating  the  motion  buffer.  This  parameter  is  called  the  Max.  Mass  Lost  During  Regrid,  and 
legal  values  are  0.0  to  0.5. 

Specifically,  the  spatial  regridding  algorithm  operates  as  follows: 

0.  Set  Total  Mass  Lost  Thus  Far  to  zero. 

1 .  Find  which  of  the  outermost  two  rows  and  two  columns  has  the  smallest  amount  of  mass. 

2.  If  this  amount  of  mass  plus  the  Total  Mass  Lost  Thus  Far  is  less  than  the  Max.  Mass  Lost 
During  Regrid  parameter,  remove  the  row  or  column  in  question  from  the  Nodestar  grid,  add  the 
mass  in  that  row  or  column  to  the  Total  Mass  Lost  Thus  Far,  and  go  to  step  1. 

3.  The  new  Nodestar  window  is  centered  in  the  remaining  area,  and  it  is  (1.0  +  Single  Step 
Motion  Threshold)  as  large  as  the  area  remaining  after  step  2  finishes.  This  enlargement  is 
performed  to  handle  the  next  motion  step.  The  new  window  will  have  the  same  number  of  cells  as  the 
original  window  and  the  mass  in  the  original  grid  is  distributed  to  the  new  grid  in  proportion  to  the 
area  of  overlap  between  each  of  the  corresponding  cells. 

4.2.2  Velocity  Regridding 

The  velocity  grid  is  in  the  standard  (r,0)  planar  coordinate  system.  The  coordinate  r  represents 
the  speed  of  the  target  as  distance  from  the  center  and  0  its  heading  in  degrees.  The  Operator 
Machine  Interface  analog  of  the  spatial  window  is  the  velocity  wheel  shown  in  Figs.  3  and  4.  The 
extent  of  the  wheel  is  set  at  the  beginning  of  the  problem  by  specifying  a  minimum  speed  rj  and  a 
maximum  speed  r2.  The  velocity  wheel  is  the  annulus  of  (r,0)  points  that  lie  between  ri  and  r2.  As 
with  the  spatial  grid,  the  number  of  velocity  cells  stays  constant  over  tiirie  for  a  given  track  (target). 

Regridding  velocity  cells  is  a  more  delicate  problem  than  regridding  the  spatial  cells.  In  the  case 
of  the  spatial  grid,  target’s  position  will  change  slowly,  so  that  Nodestar  can  move  the  Nodestar 
window,  keeping  the  high  probability  area  in  the  center.  If  the  target  starts  to  move  out  of  the  current 
window  area,  it  will  do  so  slowly  enough  that  Nodestar  will  be  able  to  shift  the  window  so  that  it 
continues  to  contain  the  target. 

By  contrast,  the  target  can  change  its  velocity  quickly.  If  Nodestar  temporarily  abandoned  some 
set  of  velocities,  then  it  would  be  quite  possible  for  the  target  to  select  one  of  these  when  it  changes 
velocity.  In  this  case.  Nodestar  might  lose  the  track.  Nevertheless,  it  is  important  to  regrid  the  set  of 
velocities.  The  way  that  Nodestar  does  this  is  to  split  just  one  velocity  cell  at  a  time,  or  combine  just 
two  cells.  In  the  spatial  regridding,  all  cells  are  resized  at  once. 

Each  velocity  cell  is  created  from  another  by  a  splitting  operation;  one  can  think  of  the  set  of 
velocities  as  a  binary  tree  with  the  root  being  the  original  velocity  cell  (containing  all  conceivable 
velocities).  There  are  two  operations: 

1.  Nodestar  can  split  a  velocity  cell  into  two  cells.  These  cells  are  called  “child”  cells,  and 
initially  each  of  them  is  “childless.”  Furthermore,  they  make  up  a  pair  of  “siblings.” 

2.  Nodestar  can  combine  two  childless  siblings,  making  the  (common)  parent  childless. 
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Suppose  that  there  are  k  velocity  cells.  Nodestar  starts  with  the  original  velocity  cell  (i.e.,  the 
annulus)  and  performs  the  splitting  operation  k  -  1  times.  At  each  time  period,  Nodestar  calculates 
the  marginal  distribution  on  velocity.  This  provides  a  probability  for  each  cell  in  the  wheel.  In 
addition  to  the  fixed  number  k  of  velocity  cells  there  is  a  threshold  probability.  After  each 
computation  of  a  posterior  target  state  distribution,  Nodestar  reviews  the  set  of  velocities  and 
combines  any  childless  pair  of  siblings  whose  total  probability  does  not  exceed  the  threshold. 
Suppose  that  Nodestar  combines  n  pairs.  It  then  splits  the  n  most  likely  cells.  Cells  can  be  split 
radially  or  angularly.  In  Fig.  3,  Nodestar  starts  with  a  velocity  wheel  and  combines  two  pairs  of 
childless  siblings  in  step  1  and  step  2.  In  Fig.  4,  Nodestar  splits  two  other  velocity  cells. 

4.3  Motion  Model 

In  Nodestar,  the  target  state  process  {X(t);  t  >  0}  is  Markovian  with  transition  function  p  defined 
by 


Vi)  =Pr  {X(t,)  =  x,|x(t,_,)  =  x,_,}  for  e  X,  and  /  =  1,2. 


(23) 


so  that 


;?(x„...,x„)  =  ,Xi.i)Po(Xo)dXo  . 


(24) 


X  <=i 


Since  the  state  space  is  discrete,  the  integral  in  Eq.  (24)  must  be  interpreted  as  a  summation. 


Initial  Wheel 


Step  1 

Fig.  3  —  Combining  velocity  cells 


Step  2 
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Initial  Wheel  Step  1 

Fig.  4  —  Splitting  velocity  cells 


Step  2 


4.3.1  Motion  Update 

Nodestar  uses  the  recursive  method  of  computing  the  posterior  distribution  on  target  state. 
Suppose  that  we  have  computed  q(tn-lr),  the  posterior  distribution  on  target  state  at  time  tn-i  .  The 
recursion  for  computing  the  posterior  distribution  on  target  state  at  time  tn  divides  into  two  steps,  a 
motion  update  and  an  information  update.  We  apply  the  Markov  chain  motion  model  to  compute  q' 
(tfi,-),  the  motion-updated  distribution  at  .  If  the  state  space  has  K  cells,  then  q(tfi.],-)  becomes  a 
column  vector  with  K  components  and  becomes  a  Khy  K  matrix  of  transition  probabilities  where 
Pn(i>j)  is  the  probability  that  the  target  will  transition  to  cell  i  at  time  r„  given  it  is  in  celly  at  time  tn-i- 
In  this  case  we  could,  in  theory,  compute  q'(tn,-)  by  matrix  multiplication  as  follows: 


—  Pnqifn-U')-  (25) 

However,  if  K  is  large,  this  computation  becomes  very  unwieldy.  Nodestar  uses  an  alternative 
method  to  perform  this  computation  efficiently. 

The  Nodestar  tracker  uses  an  underlying  Markov  chain  target  motion  model  called  a  generalized 
random  tour.  Suppose  the  target  starts  at  state  x  =  (z,v,a)  where  z  is  position,  v  is  velocity,  and  a  is 
attribute  (e.g.,  target  class).  It  moves  with  velocity  v  for  a  random  time  x  having  an  exponential 
distribution  with  mean  1/X(x).  At  this  point,  the  target  "chooses"  a  new  velocity  from  a  distribution 
that  is  a  function  of  x+vx=  (z+vx,  v,  a)  and  t.  The  target  continues  moving  with  this  new  velocity  for 
an  independent,  exponentially  distributed  time,  and  then  the  process  repeats.  Nodestar  uses  a 
discrete-time  version  of  this  process  in  which  a  new  velocity  is  drawn  at  time  tn-j  with  probability, 
l-expI-Afxjf'r^-r^.;)],  and  the  velocity  is  constant  during  (r«-7,  ?«]. 

To  describe  the  Nodestar  computation,  we  introduce  the  following  notation.  Let 


^r„(v  \z,  O')  =  Pr{  velocity  =  v  at  time  I  tgt  pos  =  z  &  attr  =  a  at  time  f„}, 
h,^{v\z,a^=  Prjnew  velocity  =  v  I  "choice"  is  made,  tgt  pos  =  z  &  attr  =  a  at  time 
r,^{y\z,a^  =  Pr{ velocity  =  v  during  I  tgt  pos  =  z  &  attr  =  a  at  time 

To  obtain  q'(tn,-},  we  let  An  =  /«  -  tn-h  and  for  each  (z,v,a)  and  compute 
r,  (v\  z,  a)  =  (y\z,a)  +  (l-  (v\z,a) 


(26) 
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? '  ,  (z  +  5(^„ ,  v)),  V,  a  J  =  (v  |z,a)2  ^(4-p(zy  .a)) 


(27) 


Note  that 


g,.(v|z,a)  = 


qitn^{z,v,a)) 

'^q(tAzy,a)) 


for  {z,v,a)  e  X 


(28) 


is  computed  after  the  information  update  at  time  The  function  6  determines  whether  there  is  a 
translation  of  the  probability  in  cell  (z,v,a)  to  another  cell.  Suppose  that  v  =  (vi,v2)  and  that  the 
spatial  grid  has  cells  of  width  di  in  the  ith  direction.  For  each  velocity  v,  the  program  recursively 
computes  an  offset  ci(v,tn)  for  the  ith  dimension.  At  the  beginning,  we  set  ci(v.O)  =  0.  After  the  nth 
motion  update,  ci(v,tn)  rs  incremented  by  ^n^i-  If  ci(v,tfi)  >  di,  then  the  ith  component  of  b(v,tn) 
equal  to  kdt,  where  k  is  the  largest  integer  such  that  kdi  <  ci(v,tn)  ^**<1  ci(v,tn)  is  decremented  by  kdi. 
This  process  represents  an  approximate  numerical  integration  of  the  velocity  process  to  compute  the 
spatial  displacement  of  the  target. 

Spherical  Coordinates.  Appendix  A  discusses  how  Nodestar  adapts  this  motion  model  to 
spherical  coordinates. 

Information  Updating.  Nodestar  performs  its  information  update  through  the  use  of  likelihood 
functions.  Having  described  the  state  space  and  motion  model  used  in  Nodestar,  it  remains 
only  to  describe  the  likelihood  functions  that  have  been  implemented  in  Nodestar  to  complete 
the  description  of  the  Nodestar  Bayesian  filter.  The  next  section  describes  the  likelihood 
functions  that  have  been  implemented  for  the  Spotlight  version  of  Nodestar. 


5.  NODESTAR  LIKELIHOOD  FUNCTIONS 

In  this  section,  we  describe  the  likelihood  functions  that  have  been  implemented  in  the  spotlight 
version  of  Nodestar. 


5.1  Bearing  Likelihood 

Nodestar  uses  the  following  method  for  computing  likelihood  functions  from  bearing 
observations  with  time-correlated  measurement  errors. 

5.1.1  Model  for  Errors  in  Bearing  Observations 

Let  Qfi  and  0^-1  be  two  bearing  observations  from  a  sensor  at  times  and  tn-h  respectively.  Let 
Bn  and  be  the  actual  bearings  to  the  target  at  times  and  tn-h  respectively  (see  Fig.  5.)  Note  the 
bearings  are  calculated  in  spherical  coordinates.  Let 

e„  =  Qn-Bn  and  En-l  =  %-l-Bn-l 

be  the  errors  in  the  two  bearing  observations.  We  assume  that  the  errors  (£«,£«- 1)  have  a  bivariate 
normal  distribution.  Specifically 
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where 


<  P<^n(^n-l 

p(Jn(7n  -  1  ^fi-i  J 


In  other  words,  and  E„-i  are  normal  random  variables  with  mean  0,  variance  of  and  (T^_,  ,  and 
correlation  p. 

For  calculating  bearing  likelihood  functions  in  Nodestar,  we  make  the  following  assumptions: 

1 .  Suppose  that  we  have  a  sequence  of  bearing  observations  0i,02,...,  0^  at  times  <  t2  <  -  < 
with  errors  ei,  ...  ,  e„.  Then  the  conditional  distribution  of  e„  given  ei,  ...  ,  Zfi-\  depends  only  on 
En-i,  i.e.,  the  bearing  errors  form  a  Markov  process. 

2.  The  correlation,  p,j  between  e„  and  depends  only  on  the  time  difference  -  t«-i.  Thus, 
there  is  a  correlation  function  c  such  that  p«  =  c(t„  -  t«-i). 


Sensor  at  tn-i 


Fig.  5  —  Observed  and  actual  bearing  to  target 


Let 


Lfi(£„  I  e„-i, ... ,  El)  =Pr{  obtaining  bearing  error  £„  I  errors  £„.i, ... ,  Ej }. 

By  assumption  1  we  have 

Lfi(£„  I  Zn-X,  ...  ,  ei)  =  Lb(£„  I  Zn-x)-  (29) 

Since  (£„  I  £„.i)^  is  bivariate  normal  with  correlation  c{tn-tn-x),  it  follows  from  Morrison 
(1976,  page  92)  that 


(£„l£„.i)  ~  iV(p£„.ia„/cy„.i,(l-p2)(j„2) 
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where  p  =  c(l„  -  tn-l)-  Let  <))  be  the  density  function  of  a  unit  normal  random  variable,  then 


VeJe„.J  =  (l) 


(30) 


5.1.2  Application  to  Nodestar 

Let  Lb(8i)  be  the  likelihood  function  for  the  bearing  error  ei.  Then 

Lb(£i)  ~  N  (0,(51^). 

By  using  Eq.  (29),  we  may  write  the  joint  probability  Lj5(e„  I  En-i,  ...  ,  ei)  of  obtaining  errors 
£l, ... ,  tn  as 

Z.,(£„...,£.)  =  L.(£,]lit(£,.„|e,).  (31) 

5.1.3  Joint  Likelihood  Function 

Suppose  that  xj,  i  =!,...,«  is  the  target  state  at  times  t,,  i  =l,...,n.  For  Nodestar  we  assume  that  the 
target  state  is  six-dimensional  and  the  motion  process  is  Markovian.  Let  p  be  the  transition  density 
for  that  process  and  pQ  be  the  distribution  on  target  state  at  time  tQ. 

Suppose  line  of  bearing  detections  are  made  at  times  ti,t2<—Jn  ^^d  that  0i,... ,  Qn  are  the  observed 
bearings.  Let  Bi(x)  be  the  true  bearing  from  the  sensor  to  the  target  at  time  ti  given  the  target  state  is 
X.  By  Eq.  (31)  and  the  Markov  assumption  on  the  target  motion  process,  one  can  write  the  joint 
probability  (density)  of  the  target  being  in  states  xq,  xi,... ,  x„  and  having  bearing  errors  ei  =  0j-Bj(x/) 
at  times  t\, ... ,  tn  as 


/*r{(Xo,Xj,...,Xjj)  and  (£j,...,e„)} 

=  PQ{Xo)f\Pi{x^\x,_x)Lg{e^)Lg{e^  £,._,) 


(32) 


<=i 


=  Po(^o)A  x^)Lgie^)UPi  iXi\x,_^)Lg(e^  £,_,). 


1=2 


An  important  feature  of  Eq.  (32)  is  that  it  allows  us  to  recursively  update  for  the  line  of  bearing 
observations. 


5.1.4  Recursive  Updating 

Let  q-{tn,-)  be  the  probability  distribution  on  target  state  obtained  from  information  updating 
through  time  tn-\  and  motion  updating  through  time  Note  that  the  bearing  error  depends  only  on 
target  position  and  not  on  the  other  components  of  target  state. 

Let  X  be  the  target’s  state  at  time  Then  x  =  {z,v,a)  where  z  is  the  target's  position  (e.g.,  lat,  Ion, 
depth),  V  is  the  two-dimensional  velocity  vector,  and  a  is  an  attribute.  For  this  calculation,  we  assume 
that  the  target  does  not  change  velocity  during  the  interval  tn-lJn-  Note:  if  Xix)(tn-tn-l)  is  small  (e.g., 
X(x)  =  1/60  minutes  and  tn-tn-l  =  10  minutes),  this  is  a  reasonable  assumption.  If  tn-tn-l  is  large,  then 
c(tn-tn-l)  becomes  small  and  the  assumption  does  not  affect  the  likelihood  calculation. 


Compute 
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^n  —  ^n~  ^n-1 


e„  =  0„  -  Bnix)  =  0„  -  Bn(z,v,a) 


En-l  =  0/1-1  -  Bn-l(x-AnV)  =  0„.i  -  5„.i(z-A„v,  V,  a) 


P  —  c{An)‘ 


Note  that  5«.i,  e^,  and  e„.i  depend  on  x  and  on  the  sensor's  position  at  times  and  tn-i-  Also 
observe  that  under  our  constant  velocity  assumption,  z-A«v  is  the  target's  position  at  time  r„.i.  Now 
we  can  compute 


l(x)  = 


>>  =  <!> 


-pg«-lgn/g„-l' 


forj^e  X 


(33) 


where  X  is  the  set  of  possible  target  states.  Finally  we  update  for  the  bearing  observation  by 
multiplying  I  pointwise  by  q'{tn,-)  to  obtain 


l(x)q~(tn,x)  forxeX. 


The  standard  deviations  0^  and  o«_i  will  be  supplied  as  part  of  the  bearing  detection  report.  For 
the  function  c,  we  use 


c(5)  =  for  5  >  0 

where  T  is  the  integration  time  of  the  sensor. 

5.1.5  Mirror  Bearings 

Mirror  bearings  are  easily  handled  in  Nodestar.  If  there  is  an  unresolved  mirror  bearing  at  time 
tn,  that  means  the  observation  at  time  tn  is  either  or  9^  where  0^  is  the  mirror  bearing  to  0«.  Let 

/  and  be  the  likelihood  functions  corresponding  to  the  observations  0„  and  9^  ,  respectively. 

tn 

Then  the  likelihood  function  /'  for  the  observation  either  0^  or  0^  is  given  by 

l'{x)  =  l{x)  +  l^{x)  ioxxGX.  (34) 

5.2  Three-Dimensional  Bearing  Likelihood 

Normally  when  incorporating  a  towed  array  bearing  into  a  Nodestar  map,  two-dimensional 
geometry  is  used  to  compare  the  observed  bearing  with  the  bearing  implied  by  a  given  target  state. 
This  is  the  case  for  the  bearing  likelihood  functions  presented  above.  When  the  target  and  array 
depths  are  small  compared  to  the  ranges  involved,  this  two-dimensional  approximation  is  reasonable. 
Two  exceptions  are  bearings  that  result  from  detection  of  a  surface  ship  near  an  array  and  bottom 
bounce  detections  of  a  submarine  by  an  array.  In  this  section  we  describe  how  one  can  account  for 
three-dimensional  bearings  in  the  Nodestar  bearing  likelihood  function. 
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5.2.1  Bottom  Bounce  Bearings 

Figure  6  shows  the  depression/elevation  (DE)  angle  resulting  from  a  bottom  bounce  path  from  a 
target  at  range  R  from  the  array.  We  make  the  idealized  assumptions  that  the  sound  rays  travel  in 
straight  lines  and  that  the  angle  of  incidence  is  equal  to  the  angle  of  reflection  at  the  ocean  bottom. 
Figure  6  shows  how  to  calculate  the  DE  angle  in  terms  of  range,  depth  of  target,  depth  of  the  receiver, 
and  ocean  depth.  The  figure  also  shows  that  the  apparent  depth  of  the  target  from  the  array  is 
d  =  2D-zrZ3,  where  D  =  ocean  depth,  Z3  =  target  depth,  and  Zr  =  receiver  depth. 

Consider  a  target  that  is  at  real  (or  apparent)  depth  Z3  and  whose  position  projected  up  to  the 
plane  of  the  array  is  (zi,Z2)  in  the  coordinate  system  shown  in  Fig.  7.  The  three-dimensional  range 
and  bearing  of  the  target  from  the  array  are  denoted  by  r  and  0.  The  two-dimensional  bearing  of  Ae 
target  in  the  (zi,  Z2)  plane  is  denoted  by  P]-.  If  the  towed  array  detects  this  target,  it  will  call  a  bearing 
6  +  e  where  e  is  the  error  in  the  observed  bearing.  If  Z3  is  small  then  6  =  pr-  In  this  case  one  can 
simply  deal  with  two-dimensional  bearings  when  computing  the  likelihood.  If  Z3  is  large,  as  it  will  be 
for  the  apparent  depth  d  resulting  from  a  bottom  bounce  detection,  then  one  must  account  for  the  DE 
angle  in  computing  the  likelihood.  We  propose  the  following  method. 


Surface 


Pig^  6  —  Depression  angle  computation  for  bottom  bounce  paths 

Establish  a  span  of  ranges  where  bottom  bounce  paths  are  possible.  Suppose  this  span  is  [/?i,/?2]- 
For  each  target  state  x,  where  x  is  a  six-dimensional  vector  (z,v,a)  (where  z  =  (zi,  Z2>  ^3))  in  the 
Nodestar  state  space,  calculate  the  two-dimensional  range  r{x)  of  the  state  from  the  towed  array.  If 
r(jc)  is  not  in  [i?i,/?2].  then  calculate  the  likelihood  as  described  earlier.  If  <  r{x)  <  R2,  then 
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compute  the  apparent  depth  d  =  2D-Zr-Z3  and  the  apparent  angle  0(j:)  from  the  relationships  in  Fig.  7, 
i.e.,  cos‘1 


.  d(x)  =  cos  ' 


Jzf+Z2+zl  J 


(35) 


One  then  compares  0(;c)  to  the  observed  angle  0  in  calculating  the  Nodestar  likelihood.  If  for 
example,  the  observation  error  8  is  normal  with  mean  0  and  standard  deviation  a,  then  the  likelihood 
for  state  x  is  equal  to 


L(x)  = 


420  [  2cr^  j 


Array 


r  =  three-dimensional  range  from  array  center  to 
s  =  two-dimensional  range  to 
0  =  three-dimensional  bearing  to 
Pr  =  two-dimensional  bearing  to 

DE  =  Depression  /  elevation 

COS0  =Z2/r,  cos  pr  =  Z\  /s,  cos  (DE)  =  sA*,  cos0=  cospr  cos  (DE) 

Fig.  7  —  Three-dimensional  angle  measured  by  an  array 
5.2.2  Bearings  to  Surface  Ships 

If  one  is  tracking  a  surface  ship  with  an  array  (fixed  or  mobile)  that  is  at  depth  D,  then  one 
calculates  0(x)  using  Eq.  (35)  with  z3  =  D.  Observe  that  when  is  small  compared  to  +Z2  .  the 
three-dimensional  bearing  is  approximately  equal  to  the  two-dimensional  bearing. 
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5.3  Detection  /  No  Detection  Likelihood 

This  section  describes  the  likelihood  function  that  is  implemented  in  Nodestar  to  represent  the 
detection/no  detection  information  from  a  passive  acoustic  sensor  such  as  a  towed  or  stationary  array. 

Consider  a  single  (moving  or  stationary)  sensor  and  let 


X  =  (  Xn) 

be  the  vector  of  target  states  at  times  ti,  i  =  1, n  and 

D  =  (Du  Dn) 

be  the  vector  of  detection/no  detection  states  at  the  sensor  at  ti,  i  =  1,  ....  n.  We  set  D,-  =  1  or  0  to 
represent  detection  or  nondetection  at  time  t,.  In  this  section,  we  develop  a  method  for  computing 
Pr{  DIx},  i.e.,  the  probability  of  observing  a  given  set  of  detections  and  nondetections  at  the  sensor 
conditioned  on  the  target's  path  x. 

5,3.1  Detection  Model 

The  model  that  we  use  for  detection  is  the  discrete  time  A,-o  process  developed  by  J.  D.  Kettelle 
and  first  studied  extensively  in  Loane  et  al.  (1964). 

For  this  model,  we  assume  that  from  the  position  of  the  sensor  and  the  target  at  time  fj,  one  can 
compute  an  instantaneous  probability  of  detection  pi  for  t  =  1,  ....  n.  This  is  typically  done  by 

calculating  the  mean  signal  excess  SE  i  in  dB  from  the  passive  sonar  equation,  e.g., 

=  SL-TL-LE-  DT, 


where 

SL  = 
TL  = 
LE  = 
AN  = 
DI  = 
SN  = 
AG  = 
DT  = 


source  level 
transmission  loss 
(AN-DI)®(SN-AG) 
ambient  noise 
directivity  index 
self  noise 
array  gain 

detection  threshold  (or  recognition  differential) 


and  ©  indicates  power  sum. 


In  the  spotlight  version  of  Nodestar,  SEi  is  computed  by  using  the  SPARS  performance 
predictions  for  the  stationary  and  mobile  arrays  of  the  Integrated  Underwater  Surveillance  System 
(lUSS).  SPARS  provides  the  LE  and  TL  terms.  Nodestar  supplies  the  DT  and  SL  terms.  DT  is 
obtained  from  an  internal  table  that  gives  DT  (also  called  recognition  differential,  RD)  as  a  function 
of  the  signal  processing  resources  that  are  on  the  beam  covering  the  geographical  location  of  the  state 
X  and  the  designated  frequency.  If  there  has  been  a  detection,  the  designated  frequency  is  the  one  on 
which  the  detection  has  taken  place.  If  there  has  been  no  detection,  then  the  frequency  is  a 
narrowband,  speed-independent  frequency  designated  for  the  class  of  submarine  given  in  the  target 
state.  Recall  that  target  class  is  the  sixth  dimension  of  the  target  state  in  Nodestar.  The  source  level 
term  5L  is  obtained  from  a  threat  database  that  provides  estimated  source  levels  as  a  function  of 
frequency  for  each  target  class. 
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We  assume  that  the  actual  signal  excess  SEi  is  given  by 


SEi  =  SEi  +  ^i 


(36) 


where  ~  N  (0,  a^).  (Note,  N(0,  a^)  denotes  a  normal  random  variable  with  mean  0  and  variance 
a^.)  A  detection  occurs  when  SEi  >  0.  Thus 

Pi  =  Pr{SEi>0}.  (37) 

In  order  to  compute  the  joint  probability  of  a  sequence  of  observations,  i.e.  D  =  (Di,  D„),  we 
must  specify  the  time  dependent  correlation  in  the  SEi’s.  That  is,  we  must  specify  the  stochastic 
process  that  yields  the  sequence  {^i,  i  =  1, ...,  n)  in  Eq.  (36).  This  is  done  as  follows:  Let  {Fq,  Y\,  ...} 
be  a  sequence  of  independent,  identically  distributed,  A(0,  o^)  random  variables.  At  each  time 
increment  ti  to  6  +  1.  there  is  the  probability 

P  =  1  -  exp  (-XA)  (38) 

of  having  a  "jump."  Let  (V,  =  the  number  of  jumps  in  [0,  6]-  Define 

Yo 

=  Yni  for  1=1,2,  ...  . 

When  a  "jump"  occurs,  an  independent  value  of  is  chosen.  In  this  case,  we  say  there  has  been 
an  acoustic  fluctuation.  These  fluctuations  are  used  to  model  the  random  fluctuations  in  SNR  that  are 
observed  for  underwater  acoustics.  Typically  X  =  1/hr  for  a  moving  sensor  and  o  is  in  the  range  of 
4  to  9  dB. 

With  this  model,  the  event  {D  i  =  1,  Z)2  =  0,  D3  =  1}  is  equivalent  to  the  event 
[SE\  >  0,  SE2  ^  0,  SE-},  >0}.  Loane  et  al.  (1964)  is  concerned  with  calculating 
Pr{Di  =  0  for  i  =  1,  ...,  n),  i.e.,  the  probability  of  no  detection  in  the  interval  [0,  tn\.  One  minus  this 
probability  is  called  cumulative  detection  probability  (cdp).  Thus,  cdp  is  the  probability  of  at  least 
one  detection  in  the  interval  [0,  6]-  Appendix  A  of  Curtis  et  al.  (1966)  gives  an  algorithm  for 
calculating  cdp  for  this  situation.  However,  this  algorithm  is  limited  to  the  case  where  Di  =  0  for 
t  =  1, ...,  n  and  is  too  computationally  intensive  for  use  in  Nodestar. 

5.3.2  Recursive  Algorithm 

For  Nodestar,  we  will  use  the  recursive  algorithm  described  below.  This  algorithm  is  based  on  the 
approximation  of  assuming  that  Di  is  conditionally  independent  of  and xi-n  forn  >  2  given  D/.i 
and  xj-i.  To  calculate  Rr{D  I  x},  we  calculate  the  four  conditional  probabilities  given  below: 

b(0,  0,  Xi,  Xi.i)  =  Pr{Di  =  0  I  D/.i  =  0,  Xi,  Xi.\ } 


b{\,  0,  Xi,  Xi.\)  =  Pr{Di  =  1  I  Di.\  =  0,  Xi,  x/.i } 


(39) 


b(0,  1,  Xi,  Xi-i)  =  Pr{Di  =  0  I  Di.\  =  1,  Xi,  Xi.\ } 


b{\,  1,  Xi,  Xi.\)  =  Pr{Di  =  1  I  Di.\  =  1,  Xi,  x/.i }. 
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To  calculate  &(0,  0,  Xi,  xi-i),  we  observe  that  given  the  position  of  the  target  at  times  tj  and  U-i 
(i.e.,  given  Xi  and  JCj+i)  and  knowing  the  position  of  the  sensor  at  these  times,  we  can  calculate  the 
instantaneous  detection  probabilities  pi  and  pi.\  at  these  times.  We  now  break  the  calculation  into  two 
cases. 

Case  1.  A  jump  occurred  between  ti-j  and  ti.  In  this  case,  detection  at  ti  is  independent  of 
detection  at  time  tj-i  so  we  obtain  simply 

Pr  {Di  =  0 1  Dm  =  0,  Xi,  Xi.\  and  a  jump}  =  1  -  pi- 
Case  2.  No  jump  occurred.  In  this  case,  =  ^M,  and 


Pr  {Di  =  0 1  Dm  =  0,  Xi,  Xi-i  and  no  jump} 


Pr  {^i.i<-SEi  and  Ul  ^  -  SEj.i} 
Pr  {^i-i  <  -  fti-i) 


=  min  {(1  -pi),  (1  -  pm)  }  /  (l-p/-l)- 

Since  P  is  the  probability  of  a  jump  occurring  in  fM  to  ti,  we  can  combine  the  two  cases  to  obtain 


\  /  \  /  xmin(l-p,),(l-p,._i) 


1-A-i 

(40) 

In  a  similar  fashion,  we  can  compute 

j/in  \  n  Ml 

(41) 

i/ni  ^  Pn  nMH 

&(0,l,X;,x,._i)-p(l-p,)  +  (l  P) 

Pi-i 

(42) 

Pi-\ 

(43) 

In  addition,  we  define 

bi(Di,xi)=  PI  if  Di  =  1, 

1  -  pi  if  Di  =  0. 

(44) 

From  Eqs.  (40)  to  (44)  it  is  evident  that  b  and  bi  depend  on  xi,  xi.\  and  the  sensor  position  at 
times  ti  and  tf^i,  only  through  pi  and  For  convenience,  we  write  an  alternate  form  of  Eqs.  (40)  to 
(44)  as  a  function  of  Df ,  pi ,  and  pi-i.  Specifically,  let 
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(40') 

(41’) 

(42') 

(43’) 


=;^1  if  Di  =  l  (44') 

\-p\  if  Di  =  0. 

From  the  above  discussion  we  can  see  that  this  likelihood  function  depends  on  both  position  and 
velocity.  As  an  example,  consider  the  situation  where  one  has  just  gained  contact ,  e.g.,  D/  =1,  = 

0.  Suppose  that  the  sensor  is  a  towed  array  and  that  there  is  a  convergence  zone  present  as  shown  in 
Fig.  8  (shaded  area).  This  likelihood  function  will  be  high  for  position  and  velocity  combinations 
that  have  the  target  just  entering  the  convergence  zone  (CZ)  and  low  for  those  just  leaving.  Consider 
for  example,  a  target  located  just  outside  the  CZ  annulus  at  ti-i  with  velocity  pointing  toward  the 
sensor  as  shown  in  Fig.  8.  This  position  velocity  combination  will  have  a  high  likelihood.  By  contrast 
the  same  position  coupled  with  a  velocity  pointing  away  from  the  CZ  will  have  a  low  likelihood. 
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Notes;  1.  Target  velocities  are  relative  to  ship  velocity 
2.  High  /  low  indicates  relative  likelihood  values 

Fig.  8  —  Position-velocity  dependence  of  likelihood  function  for  gain  of  detection 


We  can  now  calculate 


n 

Pr  {DIx}  =  bi{Di,xi)  n  b(Di,  D/-i,  x/,  x/.i).  (45) 

i=2 

Let  T  be  the  transition  function  for  the  Markov  process  representing  target  motion.  In  Nodestar 
this  process  has  a  six-dimensional  state  space.  We  use  xi  to  denote  these  target  states  at  time  fj.  Thus, 
we  may  calculate  pi  from  knowledge  of  xi  as  needed  for  the  above  calculation.  Let  qi  denote  the 
probability  distribution  on  target  state  at  time  1 .  Then 

n  n 

Pr  {D  and  x}  =  xi)  (xi)  fl  T{xi.i,xi)  fl  ^i-L  x/,  x/.i).  (46) 

i=2  i=2 

The  Markovian  nature  of  Eq.  (46)  allows  us  to  perform  updates  recursively  for  the  detection/no 
detection  information.  Suppose  that  qt  is  the  posterior  distribution  on  target  state  at  time  and  that 
we  observe  detection  information  Dj+i  at  time  f/+i.  The  update  proceeds  in  three  steps: 
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1.  Motion  Update.  Calculate 

ny,x)  for  X  G  X  (47) 

3'eX 


where  X  is  the  set  of  all  possible  target  states. 

2.  Information  Update. 

^t+l  ~  ^7+ for  X  e  X.  (48) 

3.  Renormalize  to  obtain  ^j+i-  In  Eq.  (48)  it  may  appear  that  we  need  to  know  more  than  the 
state  X  at  time  t/+i  to  compute  b.  But,  in  fact,  the  state  at  time  ti+i  contains  the  velocity  of  the  target, 
and  the  motion  model  assumes  that  the  target  holds  a  constant  velocity  during  (r,-,  r/+i).  Thus,  we 
may  compute  the  target  position  at  time  /,•  from  the  position  and  velocity  at  time  ti+i-  As  a  result,  b 
depends  only  on  £>/+i,  Di  and  x  at  time  r/+i. 

5.4  Frequency  Likelihood 

This  section  describes  the  frequency  likelihood  function  that  has  been  implemented  in  the 
Spotlight  version  of  Nodestar. 

5.4.1  Description  of  Frequency  Likelihood 

This  likelihood  function  is  calculated  when  the  source  producing  the  detected  frequency  is 
identified  as  to  whether  it  is  a  part  of  the  propulsion  system,  a  generator,  or  some  other  mechanical 
component  of  the  submarine.  An  acoustic  contact  in  this  case  consists  of  both  a  measured  frequency 
and  an  identified  mechanical  source  of  that  measured  frequency.  The  method  for  computing  the 
frequency  likelihood  function  applies  to  both  speed-dependent  and  speed-independent  sources. 

Recall  that  the  target  state  x  =  (z,v,c)  where 

z  =  (zi,  Z2,  Z3) 
zi  =  longitude 
Z2  =  latitude 
Z3  =  depth 

V  =  two-dimensional  velocity  of  the  target 
c  =  class  of  the  target. 

In  the  threat  or  target  database,  each  target  class  c  has  associated  to  it  a  set  of  speed-dependent 
and  speed-independent  sources.  Speed-dependent  source  frequencies  (and  their  harmonics)  are 
stored  as  linear  functions  of  speed.  Sound  pressure  levels  for  these  sources  are  stored  as  functions  of 
speed  and  aspect.  Speed-independent  source  frequencies  (and  their  harmonics)  are  stored  as 
expected  min  and  max  frequencies  with  sound  pressure  levels  depending  on  aspect. 

Suppose  there  are  a  total  of  N  legitimate  source  types  available  for  each  target  class.  Let  n  denote 
the  index  running  over  source  type. 


Let 
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qtix) 

L(in,y),x) 

s{x) 

a(x) 


=  Pr{Target  state  x  at  time  t} 

=  Pr{ observing  source  n  at  frequency  y  given  target  in  state  x} 
=  speed  of  target  in  state  x 
=  aspect  of  target  to  sensor  in  state  x. 


Suppose  we  receive  a  signal  at  frequency  y  at  a  sensor  identified  as  having  been  produced  by  a 
source  of  type  n.  We  assume  that  the  measurement  error  in  the  frequency  is  normally  distributed  with 
mean  0  and  standard  deviation  Of.  Let  Oa  be  the  standard  deviation  of  acoustic  fluctuations. 


For  each  target  state  x,  we  compute  the  following: 


y„(x)  =  expected  frequency  of  source  type  n  from  target  in  state  x 

A„  (x)  =  A„  (z,v,c)  =  Doppler  shift  at  the  sensor  given  a  target  of  class  c,  position  z,  and 
velocity  v  is  radiating  at  frequency  Jnix). 


L{{n,Y),x)  =  (p 


^y-7„(x)A„(x)^ 


V 


(T 


/ 


J 


(49) 


where 


9  =  density  function  for  a  normal  (0,1)  random  variable 
and  L  gives  the  likelihood  of  observing  frequency  y  when  the  actual  frequency  is  y„(x)  +  A„(x). 

Strictly  speaking,  the  term  Of  should  represent  the  measurement  error  associated  with  the  sensor. 
Consider  the  case  of  a  speed-independent  source  with  significant  variability  in  the  expected  source 
frequency.  One  would  like  to  assign  a  high  likelihood  to  measurements  of  states  that  could  have 
produced  such  data.  If  the  source  frequency  characteristically  falls  in  a  relatively  wide  band  around  a 
particular  frequency,  it  would  be  desirable  to  increase  the  error  Of  of  the  distribution.  Thus,  a 
relatively  high  likelihood  could  still  be  assigned  to  a  measured  frequency  which  differed  (by  more 
than  just  the  measurement  error)  from  the  value  stored  in  the  database. 

The  easiest  way  to  accomplish  this  would  be  to  multiply  Of  by  some  number  p  >  1  that  depended 
on  the  variability  of  the  expected  source  frequency.  The  right  way  to  do  it  would  involve  placing  an 
a  priori  distribution  on  the  frequency  in  the  database  (with  standard  deviation  representing  something 
about  the  expected  spread  in  the  source  frequency)  and  then  computing  the  standard  deviation  of  the 
posterior  distribution  of  the  source  frequency  given  the  measurement  with  error  Of  However,  this 
approach  would  entail  other,  more  complicated  modifications  to  the  above  formulas  and  would  be  of 
dubious  utility. 

5.4.2  Implementation  of  Frequency  Likelihood  in  Nodestar 

The  frequency  likelihood  in  Eq.  (49)  is  extremely  sensitive  to  target  velocity,  particularly  if  Of  is 
small,  so  it  is  necessary  to  evaluate  this  function  more  carefully  than  any  of  the  other  likelihood 
functions  in  Nodestar.  We  have  chosen  to  maximize  the  value  of  this  function  over  the  region  of 
velocity  space  being  considered  rather  than  simply  evaluating  it  at  the  center  of  the  region.  This  can 
be  done  analytically  and  without  too  much  computational  cost. 
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5.4.3  Description  of  Likelihood  Computation 

To  maximize  the  likelihood  in  a  velocity  cell  covering  {  Tq  <  r  <  ,  0q  <  6  <  0|  },  we  can  find 

the  point(s)  at  which  /(r,  0)  =  [y  -  Yn(^)  -  A„(x)]^  is  at  a  minimum. 

The  expression  for  passive  Doppler  shift  is 


A  =  7„W 


-/J(v-v4 


•  B 


where 

B  =  Unit  vector  from  z  (target  location)  to  the  sensor 
=  Sensor  velocity 

P  =  (Speed  of  sound)"  ^  . 

Note  that  all  vector  operations  assume  a  locally  flat  Earth.  For  (v  -  )*B  «  1/p,  we  can 

approximate 

A  =  p  YnW  (v  ■  )»B 


which  is  easier  to  deal  with.  For  typical  submarine  speeds,  the  error  introduced  by  this  approximation 
is  «  1  percent.  We  will  therefore  minimize  the  function: 

f{r,  0)  =  (Y-  YnW  -  P  (v  -  )*B)^  . 

5.4.4  Speed-Independent  Sources 

For  speed-independent  sources,  Yn  is  independent  of  target  state. 

/(r,  0)  =  [y-  Yn  +  P  Yn  -B  -  p  Yn  ^  cos(0  -  Bq  )f- 
=  [A  +  Br  cos(0')]^, 

where 

A  =  (Y-Yn  +  P  Yn  v^'B) 

B  =  -P  Yn 
Bq  =  direction  of  B 

0'  =  0-Bq. 

This  function’s  only  extrema  are  the  points  for  which  r  cos(0O  =  -AIB.  Note  that  for  this  set  of 
points  f(r,  0)  =  0,  so  they  specify  the  global  minimum.  There  are  no  local  minima. 
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The  shape  of  this  curve  is  such  that  it  is  impossible  to  draw  a  finite  region  in  (r,  0)  space  whose 
interior  contains  points  on  the  curve  but  does  not  intersect  the  curve  at  least  twice.  Therefore,  no 
velocity  cell  can  contain  a  minimum  without  that  minimum  also  being  on  the  perimeter  of  the  cell. 
We  may  compute  the  minimum  over  the  cell  by  computing  the  minimum  on  the  perimeter. 

To  find  the  constrained  minimum  along  the  perimeter  of  the  cell,  we  have  to  check  the  following 
points: 

A)  The  cell’s  corners: 

{ ArQ  .  00  ^’•^''0  ’  ®1  J 

B)  Constrained  minimum  ofy[ro,  0)  in  {  0q  <  6  <  0^  }: 

i)  /(rQ  ,  Bq  ),  if  {  0Q  <  Bq  <  0j  }  and  {  A  >  -Br^  } 

ii)  ArQ  ,  Bq  +  7t),  if  {  0Q  <  Bq  +  7t  <  0j  }  and  {  A  <  Br^  } 

*  iii)y(rQ  ,  cos'^f-^j), 

if  {  -1  <  ^  <  1  }  and  {  0Q  <  cos'-^f-  <  0^  }  • 

C)  Constrained  minimum  of/Cr^  ,  0)  in  {  0q  <  0  <  0^  }: 

i)  /(r^  ,  Bq  ),  if  {  0Q  <  Bq  <  0j  }  and  {  A  >  -Br-^  } 

ii)  A^\  >  Bq  +  7c),  if  {  0q  <  Bq  +  tc  <  0^  }  and  {A  <  Br^  } 

*  iii)/(rj  ,cos-^f-^)), 

if{-l<-;^  <1  }and{0Q  <cos-^r-;^)  <0j  } 

D)  Constrained  minimum  of/(r,  0q  )  in  {  Tq  <  r  <  }: 


*  /( - - ,0o),  if{ro< - - <rj 

5cos(0o  -  ^0 )  5cos(0o  -  Bq  ) 


E)  Constrained  minimum  of  j{r,  9^  )  in  {  Tq  <  r  <  }: 


*  Note  that/ evaluates  to  the  global  minimum  in  this  case,  so  we  can  stop  computation  here  if  the  above  condition 
is  met. 
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5.4.5  Speed-Dependent  Sources 

For  speed-dependent  sources,  jn  is  proportional  to  target  speed.  We  make  the  substitution 


In  =  C^r, 

where  a  is  the  Hz  per  knot  of  the  assumed  source  type.  Then 

f(r,  0)  =  [y-  Yn(x)+  P  Yn(x)  *6  -  p  Y«(x)  r  cos(0  -  Bq  )]^ 

=  [A  +  +  Cr^  cos(0')]^  , 

where 

A  =  Y 

B  =  a  (p  ‘B  -  1) 

C  =  -pa 

Bq  =  Direction  of  B 

6'  =  e-B0  • 

As  in  the  speed-independent  case,  this  function’s  only  extrema  are  the  set  of  points  for  which /(r, 
0)  =  0,  or  when  A  +  +  Cr^  cos(0')  =  0.  And,  although  this  curve  is  more  difficult  to  visualize,  it  is 

also  infinite,  and  its  only  discontinuities  tend  toward  r  =  o°.  Using  the  same  argument  as  the 
previous  case,  we  need  only  find  the  constrained  minimum  on  the  cell’s  perimeter.  The  points  we 
need  to  consider  are: 

A)  The  cell’s  corners: 

{Aro  ,  00  ),f{ro  ,  01  IKri  ,  0o  X/'*!  -  0i  )  ) 

B)  Constrained  minimum  of/(rQ  ,  0)  in  {  0q  <  0  <  0i  }: 

i)  f(rQ  ,  Bq  ),  if  {  0Q  <  Bq  <  0j  }  and  {  A  +  Br^  >  -1  } 

ii)  /(rQ  ,  Bq  +  7t),  if  {  0q  <  Bq  +  7t  <  0^  }  and  {  A  +  Br^  <  1  } 


C)  Constrained  minimum  of  f(r^  ,  0)  in  {0q  <  0  <  0i }: 
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i)  f(r-^  ,  Bq  ),  if  {  0Q  <  Bq  <  0  j  }  and  {  A  +  >  -1  }* 

ii)  /(rj  ,  Bq  +  TC),  if  {  0Q  <  Bg  +  71  <  0j  }  and  {  A  +  <  1  } 


*  iii)/(r,,cos' 


,  (A  +  Br.) 


D)  Constrained  minimum  of/(r,  0q  )  in  {  rg  <  r  <  }: 


-4ACcos{Q„-Bg)  - B  .  . .  r  i 

2Ccos(A) 


-VB'-4ACcos(0o-Be)-B^  r  i 

*  2Ccos(6„-B.) 


2Ccos(Bo-Be) 


00  ,  if<^ro<- 


2Ccos(0o  -Bg) 


E)  Constrained  minimum  of /(r,  0^  )  in  {  rg  <  r  <  }: 


a/b^ -4ACcos(0, -Bo)  - B  _  ..r  1 

•  «  ^  ---  2CcosW-B.i^’^ 


A/B^-4ACcos(a -Bo)-B  ^  .^r 


<  ''trial  <  n 


iii)/  - 


2Ccos(a  -  Bo) 


00  ,  ifi''o<- 


2Ccos(0i-Bo) 


Once  we  have  found  the  minimum  value  of/ in  the  cell,  we  convert  it  to  a  likelihood: 


Note  that/ evaluates  to  the  global  minimum  in  this  case,  so  we  can  stop  computation  here  if  the  above  condition 


is  met. 
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L  = 


-exp 


r_fi  \ 

J  rnii 


5.5  Land  Avoidance  Likelihood 

.  Nodestar  uses  two  forms  of  the  land  avoidance  likelihood  function.  The  first  form  depends  only 
on  position.  The  second  depends  on  position  and  velocity. 

5.5.1  Position  Only  Land  Avoidance  Likelihood 

Let  La  denote  the  land  likelihood  function.  The  land  likelihood  function  that  depends  only  on 
target  position  operates  in  the  following  manner.  Let 


X  =  target  state 

=  (z,v,a)  where  z  is  position,  v  is  velocity,  and  a  is  attribute 
D(x)  =  distance  of  x  from  nearest  point  of  land 
ds  =  "shallow"  distance 

dd  =  "deep"  distance  where  dj  >  dg 

a  =  avoidance  function.  This  is  a  function  such  that 


a(r) 


Jo  for  0  <  r  < 

|l  for dj  <d<oo. 


and 

a(r)  increases  monotonically  from  0  to  1  as  r  increases  from  d^  to  d^. 


Then 


LJx)  =  a(D(x))  for x^X .  (50) 

5.5.2  Position  and  Velocity  Dependent  Land  Avoidance  Likelihood 

This  form  of  the  land  avoidance  likelihood  depends  on  the  position  and  the  velocity  of  the  target 
state.  The  velocity  dependence  also  involves  the  target  motion  model.  In  particular,  it  depends  on 
x(x)  =  l/X(x),  the  mean  time  between  velocity  changes  given  the  target  is  in  state  x.  In  this  case, 

LJx)  =  a(D(x))a{D{x-i-T(x)v))  forxGZ.  (51) 

The  intuitive  meaning  of  Eq.  (50)  is  that  we  multiply  the  position-only  likelihood  by  a  factor  that 
is  equal  to  the  avoidance  function  evaluated  at  the  position  obtained  by  projecting  the  state  ahead  at 
its  velocity  v  for  a  time  equal  to  the  mean  time  between  velocity  changes  for  that  state.  If  this 
projected  position  is  close  to  land,  then  this  state  (i.e.,  combination  of  position  and  velocity)  receives  a 
low  likelihood  value.  Thus,  position-velocity  combinations  that  tend  to  put  the  target  near  land 
become  unlikely. 

5.6  Elliptical  Contact  Likelihood 

Nodestar  can  fuse  the  standard  contact  information  that  is  provided  by  many  sensors  in  the  form 
of  a  position  and  an  elliptical  uncertainty  region.  This  is  done  by  converting  the  elliptical  contact 
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information  into  a  bivariate  normal  distribution  as  follows.  The  reported  position  {Iohq,  lato)  is  the 
mean  of  the  distribution.  The  uncertainty  region  is  defined  by  an  ellipse  that  is  specified  by  the 
length  in  nautical  miles  of  the  semi-major  and  semi-minor  axes  and  the  orientation  of  the  major  axis 
measured  clockwise  from  north.  This  ellipse  is  interpreted  as  the  two-standard  deviation  ellipse  of  a 
bivariate  normal  distribution  centered  at  (lonQ,  lato).  (Note,  since  ellipses  tend  to  be  small  in  extent, 
we  use  a  local  flat  Earth  approximation  in  the  plane  tangent  to  {lono,  lato)  for  constructing  the 
bivariate  normal  distribution.)  The  ellipse  is  converted  to  the  covariance  matrix  2  of  the  bivariate 
normal  distribution  in  the  local  flat  Earth  coordinates  with  origin  at  {lono,  lato)  and  axes  parallel  to 
the  north-south  and  east-west  directions.  Let  h(-,'L)  be  the  density  function  for  the  bivariate  normal 
distribution  with  mean  (0,0)’^  and  covariance  matrix  Z.  Let 

lat(x)  =  latitude  of  the  state  x 

lon(x)  =  longitude  of  the  state  x. 

Then  the  likelihood  function  L  for  this  contact  is 


L{x)  =  h{L{x),'L)  for.xGX, 
where 


A(;c)  = 


^60cos(/ato  ){lon{x)  -  lon^  )^ 
lat{x)  -  lat^ 


(52) 


5.7  ELINT  Likelihood 


Nodestar  has  an  ELINT  likelihood  function  that  uses  the  information  about  the  type  of  radar  that 
was  detected  as  well  as  the  locating  information  in  the  detection.  Usually  the  locating  information  is 
in  the  form  of  an  elliptical  contact.  For  this  piece  of  information  we  use  the  elliptical  contact 
likelihood  given  above.  If,  in  addition,  the  type  of  radar  detected  is  identified,  we  multiply  the 
contact  likelihood  by  the  radar  type  likelihood  function,  To  use  this  function,  we  must  have  a 
database  that  lists,  for  each  radar  type,  all  the  classes  of  submarines  (or  surface  ships)  that  carry  that 
type  of  radar.  The  radar  type  likelihood  function  I  is  defined  on  the  attribute  component  a  of  the 
state  X.  Recall  that  a  is  the  target  class.  Suppose  the  radar  is  identified  as  type  k.  Then 


l(k,x)  = 


[1  if  class  a  has  a  radar  of  type  k 

[O  if  class  a  does  not  have  a  radar  of  type  k 


for  xeX , 


(53) 


where  x  =  (z,v,a).  The  effect  of  applying  the  likelihood  function  /  is  to  set  to  zero  the  probability  on 
any  state  whose  class  a  does  not  carry  radar  type  k. 

6.  MULTIPLE  TARGET  VERSION  OF  NODESTAR 

The  Nodestar  engine  provides  for  multiple  target  correlation  and  tracking.  In  this  section,  we 
describe  the  data  association  algorithms  used  and  the  method  Nodestar  uses  for  computing  multiple 
target  state  distributions. 

Recall  that  the  process  of  data  association  is  equivalent  to  partitioning  the  set  of  observations  into 
disjoint  sets  of  associated  contacts.  For  each  set.  Nodestar  creates  a  track  object.  This  is  a  C++  object, 
i.e.,  an  object  in  the  computer  science  sense.  Each  track  object  is  an  independent  copy  of  a  Nodestar 
single-target  tracker.  It  has  all  the  properties  described  in  Sections  4  and  5.  Multiple-target  Nodestar 
is  designed  so  that  these  tracks  can  be  computed  on  a  single  CPU  or  distributed  across  multiple  ones. 
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With  this  said,  the  main  question  that  remains  to  be  answered  about  multiple-target  Nodestar  is:  how 
does  it  perform  its  data  association? 

6.1  Nonlinear  Data  Association 

Nodestar  separates  the  process  of  data  association  and  state  estimation  in  the  manner  discussed  in 
Section  2.  It  is  a  single-hypothesis,  correlator-tracker  in  the  lexicon  of  that  section.  For  each  contact, 
the  probabilities  of  association  are  calculated  in  a  Bayesian  fashion  and  the  most  likely  association 
hypothesis  is  chosen.  To  calculate  the  association  probabilities.  Nodestar  uses  a  nonlinear  Bayesian 
approach  very  similar  to  the  one  described  in  Mori  et  al.  (1986).  This  is  in  contrast  to  the  data 
association  computations  that  are  performed  for  most  Kalman-based  trackers.  Those  usually  involve 
calculating  a  geometric  missed  distance  between  two  multivariate  normal  distributions  and  applying  a 
chi-squared  test.  This  approach  is  not  Bayesian  and  is  only  appropriate  for  linear  Gaussian  trackers. 

6.1.1  Probability  Model 

Suppose  that  we  obtain  an  observation  at  time  t  and  that  there  are  N  tracks  that  exist  in  Nodestar 
at  the  time  of  the  observation.  Mathematically  each  track  is  a  collection  of  observations  that  have 
been  associated  to  a  single  target.  Nodestar  uses  these  observations  to  produce  an  estimate  of  the 
target's  state. 

When  a  contact  is  obtained,  we  assume  a  prior  probability  distribution  {r„  ;  n  =1,...,  N+l}for  the 
target  that  produced  the  contact.  Specifically,  we  set 

ri  =  probability  that  contact  was  generated  by  the  ith  target  for  «  =  1 . N 

fN+l  =  probability  contact  was  generated  by  a  new  target. 

Note  that  Nodestar  does  not  distinguish  between  new  targets  and  false  targets.  Typically, 
Nodestar  operates  with 


r^+i  =  cc  (for  some  0  <  a  <  1),  and 

=(l-a)/Nfor«  =  1 . N.  (54) 

In  effect.  Nodestar  treats  the  new  target  probability  a  as  a  parameter  that  can  be  set  by  the  user. 
We  could  derive  a  from  a  model  that  estimates  the  arrival  rate  of  new  targets  and  the  generation  of 
contacts  from  these.  Of  course,  this  model  would  require  assumptions  and  inputs.  Since  we  cannot 
predict  the  situations  in  which  Nodestar  will  be  used  as  a  tracker,  we  prefer  to  specify  a  as  a  direct 
input.  Once  a  has  been  specified,  it  seems  reasonable  to  set  the  prior  probability  of  the  contact  being 
generated  by  the  target  corresponding  to  an  existing  track  equal  to  (l-oc)/N  for  each  track. 

6.1.2  Likelihood  of  a  Response 

We  now  calculate  the  likelihood  that  a  given  target  (track)  produced  the  contact.  As  an  example, 
consider  a  line  of  bearing  detection  by  a  Sound  Surveillance  Underwater  System  (SOSUS)  array. 
This  contact  has  two  components:  the  gain  of  detection  and  the  estimated  bearing  of  the  target. 
These  two  components  comprise  the  response.  Nodestar  produces  a  joint  likelihood  function  for  this 
two-component  response. 

6.1.2.1  Existing  Tracks  (Targets) 

Let  Y  =  the  response  obtained  at  time  t,  and 


Technical  Documentation  of  Nodestar _ ^ 

L(Y  \x,n)  =  Likelihood  of  obtaining  response  Y  given  that  it  was  produced 
by  target  n  in  state  x. 

Let 

q~  (t,x)  =  Probability  of  target  n  in  being  in  state  x  at  time  t. 

Then, 

g(n)={  L(Y\x,n)q~(t,x)dx  (55) 

Jx 

is  the  likelihood  that  the  response  Y is  produced  by  target  n,  for  n  =  1,...,A^. 

6.1.2.2  New  Targets 

We  now  calculate  the  likelihood  that  the  contact  was  generated  by  a  new  target.  To  do  this,  we 
need  to  establish  some  simple  probability  models  for  the  generation  of  contacts  by  new  targets. 

6.1. 2.2.1  Line  of  Bearing  Detections  —  Suppose  the  response  is  a  line  of  bearing  detection,  i.e.,  a 
gain  of  contact  and  an  estimated  bearing.  Let  the  estimated  bearing  be  6q  with  measurement  error 
that  is  normal  with  mean  0  and  standard  deviation  o  degrees.  Suppose  that  the  distribution  of  new 
targets  is  uniform  over  a  disk  of  radius  R  so  that 

q-^ft,{r,d))  =  — —  forO<  0<36OandO</-<<i?. 

180i?^ 


Suppose  that 

L(( 00, detection) I (r,  0))  =  Pr{ response  =  (0o, detection)  |  new  target  at  (r, 0)} 

=  -^exp(~^  ---^)  for  0  <  0  <  360. 

Then 


6.1.2.2.2  Elliptical  Contacts  —  If  the  reported  contact  is  an  elliptical  position  estimate  and  the 
sensor  detection  probabilities  are  not  modeled  in  Nodestar,  then  we  compute  g(N+l)  as  follows.  We 
interpret  the  elliptical  contact  as  corresponding  to  a  position  estimate  with  a  bivariate  normal  error 
distribution  with  mean  0  and  covariance  matrix  S,  where  Z  is  derived  from  the  ellipse.  Let  /  denote 
the  density  function  for  this  distribution.  Suppose  that  the  contact  is  in  the  form  of  an  ellipse  with 
center  point  y  and  associated  covariance  matrix  Z.  Then 
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L((>',E)|;c)  =  P/-{elliptical  contact(}',i;)|  new  target  at  x} 

=  f{y-x,'L)foxx&R^. 

We  shall  assume  that  the  location  of  a  new  target  is  uniform  over  a  region  W  of  area  A  so  that 

qjt,x)  =  ^  forxG  W. 


Then 


g(N  +  l)  =  j^f(y-x,'L)dx 


~  —  for  Z  "  much  smaller"  than  A. 
A 


(57) 


6.1. 2.2.3  Association  Probabilities  —  Having  calculated  the  likelihood  of  the  response  Y  given  it 
was  produced  by  target  n,  for  n  =  1,...,  A^+1,  we  calculate  the  posterior  probabilities  of  the  contact 
having  been  generated  by  the  nth  target  as  follows: 


g(nK 


for  n  =  1,...,  A^+1. 


(58) 


Nodestar  associates  the  response  to  the  target  having  the  highest  value  of  .  In  the  Spotlight 
version  of  Nodestar,  users  can  change  or  perform  associations  manually  if  they  desire. 

6.1. 2.2.4  Creating  and  Deleting  Tracks  —  If  the  highest  value  of  occurs  for  n  =  N+\,  then  a 
new  track  (object)  is  created.  Nodestar  also  contains  logic  for  deleting  tracks.  The  user  can  set  a  time 
threshold  H.  If  no  response  has  been  associated  with  a  track  for  the  last  H  hours,  then  the  track  will 
be  automatically  deleted.  This  logic  can  be  extended  to  other  criteria.  For  example.  Nodestar  could 
be  set  up  so  that  if  the  area  of  the  86  percent  containment  region  (i.e.,  the  smallest  set  of  spatial  cells 
containing  86  percent  of  the  track  probability)  for  a  track  becomes  larger  than  a  designated  value,  the 
track  is  deleted.  The  user  can  also  delete  tracks  manually. 

6.2  Probabilistic  Data  Association 

The  Nodestar  methodology  allows  for  a  nonlinear  version  of  probabilistic  data  association  (PDA) 
that  is  more  accurate  and  natural  than  the  one  given  by  Bar-Shalom  and  Fortman  (1988).  The  Bar- 
Shalom  version  of  PDA  is  applied  to  Gaussian  distributions  and  produces  Gaussian  distributions.  As 
a  result,  it  can  not  reflect  the  multimodal  distributions  that  should  be  a  natural  result  of  PDA.  Note 
that  the  methodology  described  below  has  not  yet  been  implemented  in  Nodestar. 

Suppose  that  we  have  a  sensor  response  that  can  be  associated  with  two  tracks  m  and  n.  By  this, 
we  mean  that  sum  r„  +  r„  is  "close"  to  1,  i.e.,  the  association  probabilities  for  all  other  tracks  are 
negligible.  Let  be  the  likelihood  function  given  the  response  is  associated  to  target  n  and  the 

likelihood  function  given  the  response  is  associated  to  target  m.  For  convenience,  let  us  assume  that 
r  and  f  have  been  normalized  so  that  they  sum  to  1.  Define  modified  likelihood  functions 

L  and  L„  as  follows. 

n  m 
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Ln(x)  =  l-f^  +  f^ 


Lm(x)  =  l-r^  +  f^ 


\K(x')q'n(t>x')dx' 

X 

LJx) 

lLJx')q~Jt,x')dx' 


for  xeX. 


(59) 


We  apply  L„  andL^  to  tracks  n  and  m,  respectively.  The  result  of  applying  the  likelihood 
function  Ln  in  Eq.  (59)  to  the  distribution  q~(t,-)  is  a  posterior  distribution  on  track  n  that  is  a 
mixture  of  the  posterior  distribution  that  would  have  been  produced  if  the  response  had  been 
assigned  to  track  n  and  the  one  that  would  have  been  produced  if  the  response  had  not  been  assigned 
to  «.  A  similar  result  holds  for  the  posterior  on  track  m.  Note  the  above  process  is  easily  extended  to 
distributing  the  association  over  three  or  more  tracks. 

7.  EXTENSIONS  OF  NODESTAR 


The  Nodestar  engine  and  most  of  the  likelihood  functions  described  above  have  been  developed 
for  the  Spotlight  version  of  Nodestar.  In  this  chapter  we  discuss  some  extensions  that  have  been  or 
can  be  made  to  this  version  of  Nodestar.  These  extensions  provide  examples  of  the  flexibility  and 
extensibility  of  Nodestar.  They  are  by  no  means  the  only  ways  that  Nodestar  can  be  extended. 

The  general  process  of  extending  Nodestar  to  include  a  new  type  of  information  is  quite  simple 
to  outline.  For  any  observation  or  response  representing  a  new  type  of  information  or  sensor,  one 
must  construct  a  likelihood  function  that  computes  the  probability  of  receiving  that  response  as  a 
function  of  target  state.  This  usually  requires  two  things:  first,  a  good  physical  model  of  the  sensor  or 
response  process  and  second,  the  target  state  space  must  contain  the  right  information  for  the 
calculation,  e.g.,  velocity,  if  one  wishes  to  account  for  Doppler  in  a  frequency  reading.  If  the  state 
space  contains  the  necessary  information,  then  one  has  only  to  add  a  new  module  for  calculating  the 
likelihood  function.  Having  done  this,  the  process  of  association  proceeds  as  described  in  Section  6. 
Once  the  data  is  associated  to  a  track,  the  information  update  proceeds  by  using  Eq.  (13). 

This  process  is  very  significant.  Once  the  likelihood  function  has  been  calculated,  the  data  fusion 
process  (association  and  state  estimation)  is  performed  automatically  by  the  Nodestar  engine.  This 
makes  the  Nodestar  engine  extremely  powerful.  It  effectively  reduces  the  problem  of  fusing  the 
information  from  a  new  type  of  sensor  to  a  question  of  modeling  the  physics  of  the  sensor,  a  process 
that  is  well  understood.  Use  of  this  procedure  removes  the  mystery  and  the  need  for  ad  hoc 
procedures  from  data  fusion.  It  provides  a  general  approach  and  framework  for  performing  data 
fusion. 


7.1  Surface  Ship  Tracking 

The  problem  of  surface  ship  tracking  readily  lends  itself  to  the  Nodestar  approach.  There  are 
many  types  of  sensors  whose  information  must  be  fused  to  obtain  tracks  and  identities  for  surface 
ships.  Each  of  these  types  of  information  can  be  fused  by  using  the  likelihood  approach  described  in 
Section  6.  We  list  a  number  of  typical  sensors  and  discuss  how  information  from  them  can  be  fused 
by  Nodestar. 

State  Space.  A  reasonable  state  space  for  surface  ship  tracking  is  five-dimensional,  i.e.,  position 
(2D),  velocity  (2D),  ship  class  (or  name). 

Radar.  A  typical  radar  report  is  a  position  report  with  an  elliptical  uncertainty.  This  is  a  special 
case  of  the  elliptical  contact  likelihood  already  described. 
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Acoustic.  Acoustic  detections  may  come  from  fixed  or  mobile  arrays  of  passive  sensors.  This 
information  can  be  modeled  using  the  bearing  likelihood  functions.  If  there  is  a  ship  database  that 
provides  frequency  and  source  level  information  for  the  acoustic  signals  generated  by  the  ships,  the 
frequency  likelihood  and  the  detection/no  detection  likelihood  functions  can  be  used. 

ELINT.  If  information  is  available  by  ship  class  or  even  name,  then  the  ELINT  likelihood 
function  described  in  Section  5.7  can  be  extended  to  include  not  only  radar  type  information  but 
also  pulse  repetition  interval  and  scan  frequency. 

Visual.  Visual  detections  provide  not  only  locating  information  but  also  information  about  the 
class  and  even  name  of  the  ship  detected.  This  information  is  incorporated  into  the  distribution  on 
ship  class  or  name  represented  by  the  fifth  dimension  in  the  state  space.  If  the  information  is  good 
enough,  this  distribution  will  have  probability  one  on  a  single  class  or  ship  name  for  the  target 
detected. 

Land  Avoidance.  Land  avoidance  is  also  important  for  surface  ships.  The  land  avoidance 
likelihood  function  applies  to  this  problem  also. 

7.2  Shallow  Water  Active  Acoustic  Tracking 

In  this  section,  we  describe  an  extension  of  Nodestar  that  has  been  made  to  tackle  the  problem  of 
multistatic  active  tracking  of  submarines  in  shallow  water.  This  work  has  been  performed  as  part  of  a 
Phase  I  Small  Business  Innovative  Research  (SBIR)  project  from  Naval  Sea  Systems  Command 
(NAVSEA)  and  continues  under  Phase  II  funding.  This  version  of  Nodestar  is  called  SWATS 
(Shallow  Water  Active  Tracking  System). 

A  typical  scenario  for  the  use  of  the  SWATS  tracker  might  be  as  follows.  We  are  concerned 
about  tracking  diesel  submarines  in  shallow  water  near  land.  A  radar  contact  on  the  submarine's 
periscope  is  obtained  but  lost  as  the  submarine  dives.  Three  or  more  vertical  arrays  of  hydrophones 
are  dropped  in  the  vicinity  of  the  contact  along  with  an  active  source  pinging  at  200  to  300  Hz.  We 
wish  to  process  the  acoustic  returns  at  the  arrays  to  track  and  localize  the  submarine. 

The  shallow  water  environment  is  a  complex  and  frustrating  one  for  acoustic  detection. 
Understanding  shallow  water  detection  requires  an  understanding  of  acoustic  propagation  in  shallow 
water  environments.  This  section  presents  a  brief  discussion  of  normal  mode  acoustics  and  describes 
how  it  is  used  in  the  SWATS  tracker  to  develop  a  likelihood  function  appropriate  for  this  problem.  A 
key  feature  of  this  likelihood  function  is  that  it  uses  the  acoustic  time  series  received  at  the 
hydrophones  as  its  input.  This  provides  an  example  where  the  detection  and  tracking  function  are 
combined  in  the  tracker. 

This  section  is  split  into  two  subsections.  The  first  provides  background  information  about 
normal  mode  propagation  theory  and  how  normal  modes  are  computed.  Additional  information 
about  this  theory  can  be  obtained  from  Jensen  et  al.  (1994)  and  Tolstoy  (1993).  The  second 
describes  how  normal  modes  are  used  to  develop  the  likelihood  function  for  this  situation. 

7.2.1  Normal  Mode  Propagation  Theory 

For  shallow  water,  acoustic  propagation  is  very  complex.  As  shown  in  Fig.  9,  at  ranges  beyond 
several  depth-lengths,  ray  propagation  theories  quickly  become  overwhelmed  by  the  multiple  surface 
reflections  and  bottom  bounces.  With  increased  range,  ray  theory  quickly  becomes  inaccurate. 
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Fig.  9  —  Ray  propagation  of  sound  from  a  sound  source  in  shallow  water 


The  solution  is  to  model  the  entire  water  column  as  a  shallow  waveguide,  and  to  find  the  normal 
modes  that  satisfy  the  environmental  conditions.  Each  of  the  normal  modes  travels  in  the  waveguide 
at  differing  speeds  from  the  source  and  is  recombined  at  the  sensor  to  get  the  total  pressure. 
Typically,  higher  order  normal  modes  attenuate  quickly  with  distance  so  that  with  increased  range, 
the  first  few  normal  modes  suffice  to  describe  the  sound  propagation.  This  is  shown  in  Fig.  10. 


Fig.  10  —  Propagation  of  the  first  three  normal  modes  from  a  target 


Normal  modes  are  calculated  numerically  for  a  given  environment  and  sound  frequency. 
Nodestar  uses  the  KRAKEN  normal  mode  program  (Porter  (1991)  and  Porter  and  Reiss  (1984)) 
which  is  used  by  the  NRL  Acoustics  Branch  and  has  been  successfully  tested  against  experimental 
data. 

The  normal  mode  program  KRAKEN  uses  the  sound  speed  profile,  the  attenuation,  and  the 
density  of  both  the  ocean  water  and  bottom  to  numerically  calculate  the  modes.  The  result  is  a 
numerical  evaluation  of  each  mode  over  the  entire  depth  range  z,  and  calculation  of  the 

complex  modal  wavenumber  k^. 

Using  the  modes  calculated  by  KRAKEN,  we  compute  the  sound  pressure  (j)  resulting  from  a 
source  of  strength  Qs  (measured  in  pressure  at  a  reference  distance,  e.g.,  1  yd)  at  frequency  co  and 
depth  Zs  by 

(s^{r,z,f)  =  -inQs'^  U^{zs)  U^(z)  (k^r)  (60) 

m 

=  ^rniZs)  U^iz)  ;  (k^r)  »  1 

m 

where  t  is  time,  r  is  the  range  from  the  source,  and  z  is  the  depth  at  which  the  pressure  is  measured. 

2 

is  the  Hankel  function  of  the  second  kind  which  is  approximated  for  the  far-field  in  the  second 
line.  The  time  dependence  is  usually  dropped  in  this  equation.  Note,  in  the  far  field  limit,  that 

Hl(kr)  kr»l. 

^  Ttkr  °  Kkr 

In  the  case  of  shallow  water,  k  is  usually  fairly  large.  For  k  =  0. 5m“' ,  the  condition  kr  »  1 
implies  that  r  »  2m  or  that  the  sensor  is  much  farther  than  2  meters  from  the  source. 
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For  a  vertical  array  of  n  hydrophones  placed  at  successive  depths  z/j,  for  /i  =  1,...,  n,  we  have  that 
the  pressure  reading  at  hydrophone  h  is 

Ph  =  ^^J  ^  Qs  UmiZs)  UmiZh)  QX^i-ikmrh)  (61) 

m 

where  h  indexes  the  hydrophone  at  depth  Zh  and  distance  r/,  from  the  point  source.  The  vector  P  is 
used  to  represent  the  complex  vector  which  has  complex  components  Pi,P2.--m  Pn-  We  have 
dropped  the  time  dependence  here. 

The  Um  are  numerically  computed  normal  modes  for  the  ocean  environment  of  interest,  and  are 
calculated  using  a  normal  mode  program  such  as  KRAKEN.  The  normal  modes  are  precomputed 
once  for  each  environment.  Each  mode  has  an  associated  complex  wavenumber  km  that  is  also 
computed  in  the  normal  mode  program.  There  are  typically  30  modes  of  interest  for  a  problem; 
higher  order  modes  get  attenuated  so  quickly  that  they  can  be  neglected.  The  mode  solutions  Um  are 
usually  constant  or  very  weak  functions  of  position. 

For  a  vertical  array,  the  hydrophone  distances  rf,  are  the  same  for  each  hydrophone,  and  the  only 
hydrophone-dependent  part  of  Eq.  (61)  is  the  hydrophone  depth  Zh-  Setting  rh  =  r  for  each 
hydrophone,  we  can  rewrite  Eq.  (61)  as 

r~2^ 

Ph  =  ^  Am  Umizh)  ,  where  Am=y^  Qs  UmiZs)  exp  {-ikmr).  (62) 

m 

When  sampling  the  acoustic  field  with  a  vertical  array  of  hydrophones,  it  is  convenient  to 
represent  the  transformation  in  terms  of  the  mode  amplitudes  A  =  (Ai,  A2,  ...)  at  the  sensor  vs  the 
pressure  field  P  at  the  sensor.  In  that  case,  we  can  write  Eq.  (62)  as  a  matrix  equation:  the  pressure 
vector  components  (Pi,  P2,---,  Pn)  are  samples  of  the  acoustic  field  from  the  signal,  while  the  (Ai,  A2, 
...)  are  referred  to  as  the  acoustic  mode  coefficients  of  the  signal  at  the  sensor.  In  vector  form,  Eq. 
(62)  becomes 

P  =  EA  (63) 

where  E  =  (Ehm)  and  Ehm  =  Umizh)-  The  normal  mode  propagation  model  generates  a  set  of 
modes  for  the  signal  source.  The  modes  A  propagate  from  the  source  according  to  the  propagation 
equation  described  in  Eq.  (62).  When  the  modes  are  sampled  at  the  sensor  a  distance  r  away  from  the 
source,  we  get  the  pressure  field  P.  Most  detection  methods  look  at  detection  in  field  space  or  in 
mode  space. 

7.2.2  SWATS  Likelihood  Function 

This  subsection  describes  the  likelihood  function  for  SWATS.  We  first  describe  the  statistical 
model  used  to  derive  the  likelihood  function  and  the  method  for  dealing  with  the  ambiguity  in  target 
signal  strength.  Then  we  present  the  SWATS  likelihood  function. 

7.2.2.1  Statistical  Model  of  Measurement  Error 

We  assume  that  the  pressure  readings  are  subject  to  a  multivariate  normal  distribution  of  error 
and  that  we  can  estimate  the  parameters  of  this  distribution.  The  measurement  vector  P  is  composed 
of  complex  pressure  readings  from  each  hydrophone  in  the  array.  The  measured  pressure  is  the  sum 
of  signal  plus  noise.  We  define  the  covariance  matrix  Z  for  P  as 
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'L  =  E\ 


ppff 


pH 


(64) 


where  is  used  to  represent  the  Hermitian  transpose  of  P,  or  the  complex  conjugate  of  the  transpose 
of  P. 


Using  this,  we  write  the  multivariate  probability  distribution  for  the  measurement  P  given  a  signal 
Pq  (see  Theorem  2.3.1  of  Anderson  (1984)): 


Pr{P\F,}  =  N(P\F,p  = 


i(P-Po)^E-'(P 


(65) 


7.22.2  Signal  Strength  Ambiguity 

In  order  to  use  Eq.  (65)  as  a  likelihood  function,  we  must  be  able  to  predict  Po(x),  the  signal  at 
the  array  given  that  the  target  state  is  x.  To  do  this,  we  have  to  predict  the  pressure  at  the  array 
resulting  from  a  Qse'*'^^‘  signal  emitted  by  the  pinger  and  reflected  off  a  target  in  state  x.  Define 
po(x)  as  the  signal  (i.e.,  complex  vector  of  pressures)  at  the  array  resulting  from  the  incident  signal 
from  the  pinger  being  reradiated  from  a  point  source  at  x  with  property  that  magnitude  of  the 
reflected  pressure  at  the  reference  distance  (e.g.,  1  yard)  is  equal  to  the  magnitude  of  the  incident 
pressure  at  x. 

Define  Q(x)  to  be  the  ratio  of  the  reflected  pressure  at  the  reference  distance  to  the  incident 
pressure  given  the  target  is  in  state  x.  Let  TS(x)  be  the  target  strength  of  the  target  given  it  is  in  state 
X.  Then 


TS(x) 

Q(x)  =  lO  20  . 


Now  we  can  write 


Po(-x)  =  2(-^)  Po(^)-  (66) 

The  signal  po(x:)  is  calculated  by  using  the  modal  decomposition  to  compute  the  signal  that 
arrives  at  the  target  Pt(x),  given  the  target  is  in  state  x  and  then  assuming  that  this  signal  is  reradiated 
from  the  target  as  though  it  were  a  point  target  with  Q(x)  =  1.  The  signal  received  at  the  array  is 
calculated  by  applying  the  modal  decomposition  to  Pt(x)  and  propagating  it  to  the  array. 

7.2.2.3  Propagation  from  Pinger  to  Target 

To  calculate  the  sound  propagation  from  the  pinger  to  the  target,  we  use  Eq.  (62): 


P,ix)  =  Q^UJZp)[/Jz,(x))exp(-ik^r,(x)] 

m  \k„r,(x)  ^  > 

where 

Zp  =  depth  of  the  pinger, 
zt(x)  =  depth  of  the  target  in  state  x 
rt(x)  =  range  of  target  in  state  x  from  pinger 
Qs  =  pinger  strength  (measured  in  pressure  at  reference  distance). 
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1.2.2A  Target  Reradiation 

The  target  reradiates  the  sound  in  a  radiation  pattern  dependent  upon  the  target’s  alignment  and 
profile.  The  value  of  Q{x)  reflects  this  pattern.  For  the  purpose  of  calculating  po(^J,  we  assume 
Q{x)  =  1. 

1. 2.2.5  Propagation  from  Target  to  Sensor 

The  sensor  is  a  vertical  array  of  hydrophones  in  the  water  placed  at  depths  z/;  for  /z  =1,...,  n. 
Using  Eq.  (62)  with  Pfix)  in  place  of  Qs,  we  calculate  ph(x),  the  pressure  at  the  hydrophone  with 
depth  zh,  as  follows: 


271 


P,  iz,  (x))f/„  (z,  )exp(-z\r,  (x))  for  /z  =  1 , . . . , 


(67) 


where 


Zh  =  depth  of  the  hydrophone  h, 

Zf(x)  =  depth  of  the  target  in  state  x 

r^(x)  =  range  from  target  in  state  x  to  the  array. 

Now  we  have 


PoW  =  (PiW’p2W>--->P«W)^- 

Using  Eq.  (62),  we  rewrite  the  probability  distribution  in  Eq.  (65)  as  follows: 


Pr{Pl(2(x),Po(x)} 


i(P  -  Q(x)f,ix))'<  2-'  (P  -  e(x)p„(j:))) . 


(68) 


Equation  (68)  represents  the  probability  of  receiving  P  given  Q(x)  and  a  theoretical  unit  strength 
signal  return  pofxj. 

We  handle  our  lack  of  knowledge  of  Q(x)  by  finding  the  target  strength  gbest  that  maximizes  Eq. 
(68)  for  a  given  P  and  po(x).  To  do  this,  we  rotate  and  scale  the  covariance  matrix  Z  to  be  diagonal. 
Since  Z  is  a  covariance  matrix,  it  is  positive  definite  and  there  is  a  nonsingular  matrix  R  that  satisfies 

RZR^  =  I  =  (RZR^)  1 


Z-1  =  R^R. 


(69) 


We  can  now  find  the  target  strength  2best  that  maximizes  Eq.  (68)  which  is  equivalently 
expressed  in  Eqs.  (70)  and  (71): 

Gbest(^)=  |e:-^(P-e(^)Po(^)fR"^P-G(^)PoH)ismaximize^^  (70) 
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!2best(4:)=|!2:--(RP-!2(^)RPo(^))  (rP-<2(^)RPo(^))  is  maximizedj.  (71) 


The  target  strength  gbestC^)  that  maximizes  Pr/PI  Q(x),  in  Eq.  (71)  is  calculated  to  be 

RPo(x) 


(72) 


1. 2.2.6  Likelihood  Function 

The  likelihood  function  in  SWATS  uses  the  maximum  likelihood  signal  strength.  The  covariance 
matrix  for  P  is  assumed  to  be  so  that  the  corresponding  coordinate  transform  matrix  R  is  ctI. 

For  a  given  target  state  x,  we  calculate  po(x)  which  is  the  theoretically  expected  signal  for  a  target 
with  unit  target  strength.  We  also  calculate  \Qmax(x)\,  which  is  the  maximum  possible  target  strength 
for  target  state  x.  If  I2best(-^)l  >  I2max(x)l,  we  scale  down  gbesK^)  so  that  it  has  a  magnitude  equal  to 
Qmaxix)-  The  likelihood  function  L  that  is  used  in  SWATS  comes  from  Eq.  (68)  and  is  calculated 
within  a  scale  factor  as 


LP 


¥ 


( 


X  =  exp 


1 


P-!2„(^)PoW 


(73) 


where  the  maximum  likelihood  target  strength  Qmix)  is  calculated  as 

Po(x)^P. 


2m  (^)  =  2best(^)  = 


Po(^) 


2best(^)N 


(74) 


2m(x)  = 


2max(^) 


2best(^) 


2best(^);  GbestWP  2™axW 


(75) 


In  SWATS,  <2max('^)  is  obtained  as  an  input. 

7.3  Incorporation  of  Subjective  Information 

Subjective  information  can  be  incorporated  into  Nodestar  in  two  ways:  in  the  prior  distribution 
on  target  state  and  by  the  use  of  likelihood  functions. 

7.3.1  Prior  Distribution 

The  prior  distribution  on  target  state  contains  information  about  the  target's  initial  location, 
velocity  distribution,  and  class.  All  of  these  can  be  adjusted  on  the  basis  of  subject  (or  objective) 
information  that  is  available  before  Nodestar  begins  to  track  the  target. 

7.3.2  Likelihood  Functions 

Subjective  information  that  is  obtained  after  a  track  is  initiated  can  be  incorporated  by  using 
subjective  likelihood  functions.  A  simple  example  will  suffice  to  illustrate  how  this  is  done. 
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Suppose  that  we  have  subjective  information  about  the  class  of  one  of  the  targets  that  we  are 
tracking.  If  we  can  quantify  the  information  in  terms  of  subjective  probabilities,  then  we  can  define  a 
likelihood  function  and  update  the  target  state  distributions  accordingly. 

As  an  example,  suppose  that  we  have  K  target  classes,  /'a^' A:  =  1, ....  A'/ ,  in  the  target  state  space. 
We  might  receive  an  intelligence  report  about  the  target  that  contains  information  about  the  target's 
class.  Suppose  that  on  the  basis  of  this  information  we  assign  the  following  subjective  probabilities  to 
the  class  of  the  target: 


Pr{a  =  ai^}  =  X,^foxk-\,...,K.  (76) 

Because  information  updating  by  the  use  of  likelihood  functions  is  insensitive  to  multiplicative 

factors  (see  Eq.  (13)),  the  A^'s  need  not  add  to  one.  In  fact,  they  can  be  relative  credences  or  weights. 

Using  Eq.  (76)  we  define  the  following  likelihood  function. 

L(x)  =  L((z,v,a))  =  Xi^  fora  =  a^  andx  =  (z,v,a)G  A(.  (77) 

Let  q{t,-),  be  the  probability  distribution  on  the  target  state  at  time  t.  The  posterior  at  time  t,  given 

the  subjective  information  on  target  class,  is  computed  by 

q(t,x)  =  —L(x)q(t,x)  forxGX 
C- 

where 

C  =  ^L(x)q(t,x}dx 

X 

as  given  by  Eq.  (13). 
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Appendix  A 

SPHERICAL  MOTION  MODEL  AND  GRID 


This  Appendix  describes  the  modifications  made  to  Nodestar  to  change  from  a  flat  Earth  grid 
model  to  a  round  Earth  grid.  Section  A1  describes  the  rationale  for  the  changes.  Section  A2 
describes  the  mathematics  involved.  Section  A3  reviews  current  practice  in  applying  motion  models 
to  flat  grids.  Section  A4  describes  the  changes  we  have  made  so  that  we  can  apply  motion  models  to 
round  grids.  Section  A5  describes  the  problems  we  cannot  avoid  in  going  to  the  round  Earth  model. 

Al.  RATIONALE 

In  the  past,  Nodestar  was  a  single-target  tracker.  We  assumed  that  we  knew  approximately  where 
a  target  was  so  that  we  could  begin  the  tracking.  Thus,  our  grids  represented  only  a  small  portion  of 
the  Earth.  In  fact,  even  in  debugging  runs,  the  grids  used  in  Nodestar  never  covered  more  than  10°, 
and  so  at  most  600  nmi.  At  the  Equator,  the  amount  of  distortion  of  the  sides  of  the  grid  caused  by 
this  is  at  most  2  percent.  For  a  submarine  moving  at  10  kt  and  within  6  h  travel  from  an  island,  the 
grid  is  2°  by  2°,  and  the  distortion  is  at  most  0.01  percent. 

Because  Nodestar  is  now  a  multiple-target  tracker  designed  for  problems  involving  large  areas  of 
the  ocean,  it  must  consider  the  Earth's  curvature.  Fortunately,  the  mathematics  of  navigation  leads  a 
simple  modification  to  the  methods  used  for  a  flat  Earth.  We  can  find  the  equations  of  motion  for  a 
ship  traveling  along  a  constant  bearing,  and  our  results  are  nearly  the  same  as  the  results  we  had  for  a 
flat  Earth. 

A2.  THE  MATHEMATICS  OF  NAVIGATION 

First,  we  give  some  definitions.  The  spherical  coordinates  for  the  surface  of  the  Earth  are  ([)  and 
e,  where  (])  is  the  colatitude,  defined  as  90°  minus  the  latitude,  and  0  is  the  longitude.  Let  us  define  (p 
to  be  the  latitude,  so  that  (p  =  90°  -  (j).  Since  the  metric  for  the  sphere  in  spherical  coordinates  is  given 
by  Eq.  (Al)  below,  the  metric  in  latitude-longitude  coordinates  is  given  by  Eq.  A2  below. 

ds^  =  +  sin^  <])  dOi^,  (Al) 

ds^  =  dtp^  +  cos^  cp  dQ^.  (A2) 

A  curve  described  by  the  motion  of  a  ship  maintaining  a  constant  bearing  with  respect  to  north  is 

called  a  rhumb  line  or  a  loxodrome.  This  is  not  a  great  circle,  although  for  short  distances  it  is  close. 

A  loxodrome  spirals  into  the  north  or  south  pole.  We  will  derive  its  equation  below. 

Assume  that  a  ship  starts  off  at  latitude  cpg  and  longitude  0q,  and  that  its  equation  of  motion  is  (p 

=  (p(0,  0  =  0(0  in  such  a  manner  that  its  bearing  with  respect  to  north  is  a  and  its  speed  is  the  constant 
s.  Then,  using  the  metric  tensor  for  the  latitude-longitude  coordinate  system,  we  can  use  the  equation 
u  v  =  |ul  |v|  cos  a,  where  h  =  (1,  0)  and  v  =  (cp',  0')  to  find  the  equation 
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tp'  =  5  cos  a. 

(A3) 

To  find  0',  we  must  use  the  fact  that  the  length  of  v  is  5  to  get 

cp'^  +  cos^  (p  0'^  =  . 

(A4) 

Substituting  Eq.  (A3)  into  Eq.  (A4),  we  get 

cos^  (p  0'^  =  s'^  sin^  a. 

(A5) 

Therefore,  we  get  the  following  for  0': 

q0  =  5  sin  a  sec  j. 

(A6) 

Thus,  the  rate  of  change  of  latitude  is  the  same  in  the  round  Earth  tracker  as  in  the  flat  Earth 
tracker,  and  the  distortion  in  the  rate  of  change  of  longitude  is  quite  simple.  Incidentally,  we  can  get 
a  parametric  solution  for  (|)  and  0  by  substituting  for  s  in  Eq.  (A6); 

qi,  =  70  tan  a  sec  7, 

(A7) 

^0  cot  a  =  70  sec  7, 

(A8) 

0cota  =lntan  I2  +  4!  + 

(A9) 

A3.  FLAT  EARTH  GRIDS 

Originally,  Nodestar  used  a  simple  model  for  moving  target  positions  within  a  given  grid  and 
motion  model.  It  assumed  that  the  probability  densities  for  a  grid  cell  are  concentrated  at  the  cell 
center  as  in  point  A  in  Fig.  Al.  If  we  have  a  certain  velocity  under  consideration,  (v^  ,  ),  and  the 

grid  cell  size  is  {dx,  dy),  then  the  target  will  move  one  cell  to  the  right  at  time  tx  =  dxivx  and  one  cell 
up  at  time  ty  =  dy/vy.  It  will  move  further  to  the  right  or  up  at  every  multiple  of  tx  or  ty,  respectively. 

Quite  simply,  we  keep  track  of  the  current  time,  t  ,  and  t  ,  and  we  move  the  probabilities  for  all 

^  y 

cells  in  lockstep  at  the  multiples  of  t  and  t  Since  these  parameters  do  not  depend  at  all  on  the 

X  y 

original  cell  indices,  cells  will  never  overlap.  The  target  traveling  from  A  to  B  below,  in  our  model, 
enters  the  cells  in  the  manner  described  below. 

A  target  starting  at  point  A  traveling  toward  point  B  was  assumed  to  go  first  to  the  right,  then  up, 
then  right  again.  Each  time  that  it  passed  a  dashed  line,  denoting  that  it  passed  the  center  of  a  grid 
box  either  horizontally  or  vertically,  it  jumped  a  box  in  the  internal  representation.  The  model  would 
do  the  same  for  any  target  traveling  parallel  to  AB. 

A4.  THE  CONVERSION  TO  ROUND  EARTH  GRIDS 

The  difference  between  round  Earth  grids  and  flat  Earth  grids  is  that  the  areas  of  different  grid 
boxes  will  be  different.  The  area  of  the  grid  box  with  latitude  going  from  (p  to  cp+Acp  and  longitude 
going  from  0  to  6+  A0  is  [sin  (cp+Atp)-  sin  (p]  x  A0.  This  last  figure  is  approximately  equal  to  coscp  x 
A(p  A0,  by  the  mean  value  theorem.  The  lengths  of  the  north  and  south  sides  of  the  grid  box  are 
related  similarly. 
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Fig.  A1  —  Flat  Earth  grid 

Fortunately,  targets  change  their  latitudes  in  this  model  at  the  same  times  as  in  the  previous  model 
by  Eq.  (A3)  above.  That  equation  is  the  same  as  the  analogous  equation  in  the  flat  model.  For  the 
spherical  model,  Nodestar  no  longer  treats  every  cell  in  the  same  way.  It  does  treat  each  row  the  same 
way.  In  other  words.  Nodestar  keeps  track  of  which  row  each  target  position  was  in  at  time  0,  as 
follows. 

Define  the  local  advancement  times  for  a  velocity  model  for  a  target  traveling  at  speed  s  and 
bearing  a  in  a  grid  box  centered  at  (cp,  0)  with  subgrid  sizes  (Acp,  A0)  as  follows: 


A(p 

s  cos  a 


(AlO) 


A0 

tg  =  - COStp. 

s  sin  a 


(All) 


At  times  (q,  21q  ,3tQ  ,  Nodestar  moves  the  target  east.  At  times  t^,  2t^,  it  moves  the  grid 
box  north,  and  then  corrects  for  the  shortening  of  lines  of  latitude.  Suppose  that  at  time  t=t^, 


o<t(p  <te 


(A12) 


As  the  target  enters  the  grid  box  centered  at  ((p  +  A(p,  0),  the  value  of  will  change  to 


^6new 


A0 


s  sin  a 


-cos({p  + A(p). 


(A13) 
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To  smooth  out  the  transition  a  bit,  let  us  average  the  fg's  by  integration.  We  integrate  the 
expression  for  with  respect  to  cp,  and  we  get 


1  Ae 

A(p  5  sin  a 


cos  tdt 


(A14a) 


A0  ,  .  ,.^.sinA(p/2 

= - X  cos((p  +  A(p/2) - ^ 

5  sin  a  A(p/2 


(A14b) 


Ae 

s  sin  a 


X  cos((p  +  A(p/2). 


(A14c) 


When  a  target  goes  up  a  box,  look  at  the  time  that  has  passed  since  it  has  moved  horizontally. 
Call  this  time  th-  If  th  is  greater  than  t0avg.  we  move  the  target  right  and  subtract  feavg  from  th-  In  any 
case,  we  multiply  th  by  renew  divided  by  reavg  and  continue.  Consider  the  diagram  in  Fig.  A2. 

For  a  target  traveling  from  A  to  D,  Nodestar  moves  it  up  one  grid  line  at  point  B.  The  correction 
moves  it  to  the  right,  and  then  it  continues  to  C  and  D.  Similarly,  a  target  moving  from  E  to  I  has  a 
correction  made  at  point  G. 


A9 


Fig.  A2  —  Motion  in  round  Earth  grid 

AS.  UNAVOIDABLE  DIFFICULTIES  IN  ROUND  MOTION  MODELS 

One  might  think  that  the  problems  we  face  in  designing  homogeneous  grids  for  the  round  Earth 
model  simply  indicate  that  we  have  not  adequately  thought  through  our  problems.  This,  however,  is 
false.  No  matter  what  we  do,  we  cannot  find  equal-area,  translation-invariant  grids  for  the  sphere. 
This  is  a  matter  of  algebraic  topology,  not  of  analysis. 

Suppose  we  could  find  such  a  grid.  Then,  we  would  have  coordinates  for  the  sphere  without  any 
singularities.  We  could  get  a  nonvanishing  continuous  vector  field  on  the  sphere.  This,  in  fact,  is 
impossible.  By  the  Poincare-Hopf  Theorem  (see  Guillemin  and  Pollack  (1974))*,  any  continuous 
vector  field  must  have  at  least  one  zero  element. 


Guillemin,  V.  and  Pollack,  A.  (1974).  Differential  Topology.  Prentice  Hall,  Englewood  Cliffs,  New  Jersey. 
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This  result  means  that  any  grid  scheme  must  break  down  at  least  one  point.  We  could  create  an 
artificial  coordinate  system  on  the  sphere  where  there  was  one  pole  in  the  middle  of  a  large  land  mass, 
but  we  would  need  to  be  careful  to  have  this  be  part  of  the  internal  implementation  only. 


Appendix  B 


A  SUMMARY  OF  THE  COMPARISON  OF  LINEAR  AND  NONLINEAR 

TRACKERS 


This  report  presents  the  results  of  a  series  of  tests  designed  to  compare  the  performance  of  the 
discrete  nonlinear  tracker,  Nodestar,  to  that  of  an  extended  Kalman  filter  in  a  variety  of  situations 
involving  the  use  of  line-of-bearing  detections  to  track  a  submarine.  The  extended  Kalman  filter 
used  in  this  comparison  is  the  Maneuvering  Target  Statistical  Tracker  (MTST).  MTST  is  used  in  a 
number  of  Navy  systems  including  the  OBU  tracker  and  NTSA's  TIMS/MRS  track  reconstruction 
system  (see  Gerr  (1984).) 

Thirty  cases  were  run.  In  each  case,  there  were  three  fixed  passive  arrays  and  a  target  track  that 
was  chosen  in  a  random  fashion.  Detection  of  the  target  by  the  arrays  was  simulated  using  an  1-s 
process  for  acoustic  fluctuations.  (The  1-s  process  is  defined  in  Loane  et  al.  (1964).)  When  a 
detection  occurred,  the  simulation  calculated  the  true  bearing  of  the  target-to-detecting  sensor  and 
then  added  a  random  error.  Each  case  lasted  for  20  hours  of  simulated  time.  For  each  case,  the  file 
of  detections  and  noisy  bearings  was  recorded.  The  same  file  was  fed  to  Nodestar  and  MTST.  Every 
15  minutes  during  the  run,  each  tracker  produced  an  estimate  of  the  target's  location  in  terms  of  a 
probability  distribution.  These  estimates  were  recorded  and  used  to  calculate  measures  of 
effectiveness  (MOEs)  that  measure  the  performance  of  the  two  trackers. 

Bl.  DESCRIPTION  OF  MOEs 

Both  Nodestar  and  MTST  are  statistical  tracers  in  the  sense  that  their  estimates  of  target  position 
are  in  terms  of  probability  distributions.  MTST's  output  is  displayed  as  ellipses.  The  ellipses 
represent  bivariate  normal  distributions  for  target  location  and  are  the  contours  of  the  86  percent 
containment  region.  Nodestar's  output  is  a  discrete  probability  distribution  on  target  state.  Each  cell 
in  the  distribution  is  color-coded  to  indicate  its  probability. 

We  used  three  MOEs: 

AOU  Size.  This  is  the  area  of  the  smallest  region  containing  86  percent  of  the  probability 
distribution  on  target  location.  For  MTST  this  is  the  area  of  the  86  percent  (2-s)  ellipse. 

Mean  Missed  Distance  (MMD).  For  this  measure  we  must  know  the  target's  actual  position. 
We  then  calculate  MMD  as  the  mean  squared  distance  of  the  target's  actual  position  from  the 
probability  distribution  on  target  location.  MMD,  then,  is  the  generalization  to  probability 
distributions  of  the  mean  missed  distance  that  one  would  calculate  if  the  tracker  gave  a  point  estimate 
for  target  position. 

Accuracy.  We  measure  the  accuracy  of  the  probability  distributions  produced  by  a  statistical 
tracker  in  terms  of  its  containment  regions.  As  an  example,  take  the  86  percent  containment  region 
for  MTST  (i.e.,  the  ellipses).  Suppose  we  plotted  the  actual  target  position  at  100  times  and 
compared  it  to  the  86  percent  region  at  those  times.  We  want  the  target's  position  to  fall  inside  the 
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region  86  percent  of  the  time  and  outside  the  remaining  14  percent.  Similarly,  if  we  looked  at  the  50 
percent  region,  we  would  want  the  target  to  be  in  50  percent  of  the  time  and  out  the  other  50  percent. 
If  this  is  true,  then  our  distribution  is  giving  us  an  accurate  representation  of  the  degree  to  which  we 
have  localized  the  target. 

Accuracy  is  a  measure  of  the  extent  to  which  a  tracker's  containment  regions  are  an  accurate 
representation  of  its  knowledge  of  the  target's  location.  Accuracy  varies  from  100  percent  to  0 
percent. 

In  making  our  comparison,  we  calculated  the  AOU  size  and  MMD  at  each  time  the  trackers  made 
an  estimate  (i.e.,  81  times  for  each  case)  and  computed  the  average  AOU  size  and  MMD  as  our 
comparison  statistics.  The  accuracy  MOE  is  computed  for  an  entire  case.  It  is  based  on  the 
distributions  produced  at  the  81  times  during  the  run. 

B2.  ROBUSTNESS  TESTING 

Both  Nodestar  and  MTST  require  the  user  to  make  probabilistic  assumptions  about  the  target's 
motion.  Since  Nodestar  uses  detection/no  detection  information  from  the  arrays  to  help  estimate  the 
target's  location,  it  requires  an  estimate  of  the  figure  of  merit  (FOM)  for  the  target's  primary  (i.e., 
most  detectable)  frequency  and  the  propagation  loss  curve  for  this  frequency.  This  acoustic 
information  is  necessary  for  the  detection/no  detection  likelihood  functions  calculated  by  Nodestar. 
We  investigated  the  robustness  of  tracker  performance  when  there  is  a  mismatch  between  the 
assumptions  used  by  the  tracker  (i.e.,  the  motion  and  acoustic  assumptions)  and  those  used  in  the 
simulation  that  generated  the  file  of  detections  fed  into  the  trackers. 

B3.  TEST  SCENARIOS 

We  used  the  following  scenarios  in  our  tests. 

Target  Motion  Models.  For  the  target  motion  in  the  simulator,  we  used  three  basic  types  of 
models:  patrolling,  constant-course-and-speed,  and  evading.  We  used  the  three  base  tracks  shown  in 
Figs.  B1  through  B3.  The  figures  also  show  the  locations  of  the  three  fixed  array  sensors  (labeled  1, 
2,  and  3)  used  to  generate  detections  in  the  simulation.  The  rectangle  in  each  figure  marks  the  end  of 
the  base  track. 

Patroller.  For  the  10  patroller  cases,  we  used  the  base  track  shown  in  Fig.  Bl.  The  speeds  for 
each  of  the  legs  are  given  in  Table  Bl.  For  each  case,  we  perturbed  the  length  of  each  leg  by  making 
a  random  draw  to  adjust  the  length  of  the  leg. 

Constant-Course-and-Speed.  For  the  10  constant-course-and-speed  cases  we  used  the  target 
track  shown  in  Fig.  B2  with  a  speed  of  12  kt. 

Evader.  Figure  B2  shows  the  base  track  for  the  evader  motion.  Table  B2  gives  the  speeds  for  each 
of  the  legs.  For  each  case,  we  perturbed  the  length  of  each  leg  by  making  a  random  draw  to  adjust 
the  length  of  the  leg. 


60 


Stoney  Corwin,  and  Hofmann 


Evader  Target 


Fig.  B3  —  Evading  target  base  track 


Table  B1  —  Patroller  Tracks 


Leg 

Course  (deg) 

Speed  (kt) 

Time  on  Leg 

1 

45 

12 

RAND(2.5) 

2 

270 

11 

RAND(2.0) 

3 

0 

11 

RAND(4.0) 

4 

90 

10 

RAND(3.0) 

5 

225 

12 

RAND(l.O) 

6 

135 

12 

RAND(1.5) 

7 

225 

12 

RAND(3.0) 

8 

280 

10 

RAND(5.5) 

*Note:  RAND(x)  indicates  a  random  draw  from  a  uniform  distribution  over  [x/2,  3x/2].  The  base 
case  has  time  on  leg  equal  to  x  hr. 
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Table  B2  —  Evader  Tracks 


Leg 

Course  (deg) 

Speed  (kt) 

Time  on  Leg 

1 

no 

10 

RAND(l.O) 

2 

200 

12 

RAND(l.O) 

3 

100 

10 

RAND(l.O) 

4 

270 

13 

RAND(1.5) 

5 

220 

10 

RAND(2:0) 

6 

170 

12 

RAND(2.5) 

7 

200 

14 

RAND(12.0) 

*Note:  RAND(x)  indicates  a  random  draw  from  a  uniform  distribution  over  [x/2,  3x/2].  The  base 
case  has  time  on  leg  equal  to  x  hr. 

Tracker  Motion  Assumptions.  Both  Nodestar  and  MTST  used  one  fixed  motion  model  for  all 
30  test  cases.  The  Nodestar  motion  model  assumed  that  the  target's  heading  was  uniform  [0,  360°] 
and  its  speed  distribution  uniform  over  [5,20  kt].  The  target  course  changes  were  modeled  as  taking 
place  at  exponentially  distributed  times  with  a  mean  time  of  1  hour  between  course  changes. 

The  motion  model  in  MTST  is  an  Integrated  Ornstein-Uhlenbeck  process.  The  choice  of  the 
parameter  values  for  this  motion  model  was  made  in  conformance  with  the  recommendations  in  Gerr 
(1984).  The  objective  was  to  match  the  motion  assumptions  used  by  Nodestar  as  closely  as  possible. 

Acoustic  Assumptions.  The  test  considered  three  acoustic  environments:  surface  duct,  direct  path 
and  bottom  bounce,  and  convergence  zone. 

Nodestar  requires  that  the  user  specify  a  propagation  loss  curve  and  an  FOM  for  the  target.  To 
test  the  robustness  of  Nodestar  to  mismatches  in  the  propagation  loss  curve  and  the  FOM,  we 
considered  the  following  cases: 

1 .  Correct  transmission  loss  curve,  correct  estimate  of  FOM. 

2.  Correct  transmission  loss  curve,  high  estimate  of  FOM. 

3.  Correct  transmission  loss  curve,  low  estimate  of  FOM. 

4.  Incorrect  transmission  loss  curve,  correct  estimate  of  FOM. 

We  ran  one  test  for  each  combination  of  assumptions  listed  above.  This  produced  three  motion 
models  x  three  propagation  loss  curves  x  three  FOM  assumptions  +  three  mismatched  propagation 
loss  curves  =  30  runs. 

Each  run  was  given  a  three-letter  descriptive  label.  Table  B3  describes  the  meaning  of  the 
three-letter  label. 
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First  Letter _ _ 

Acoustic  Environment 

C  =  Convergence 
Zone 

D  =  Direct  Path 
S  =  Surface  Duct 


Table  B3  —  Run  Labels 

Second  Letter 
Simulated  Target 
Motion 
P  =  Patroller 

C  =  Constant  Course 
and  Speed 
E  =  Evader 


Third  Letter _ 

Acoustic  Information 
Accuracy  (Mismatch) 
C  =  Correct  Acoustic 
Assumption 
H  =  High  estimate  of 
FOM 

L  =  Low  estimate  of 
FOM 

I  =  Incorrect 
Propagation 
Loss 


B4.  RESULTS  OF  COMPARISON 

Table  B4  summarizes  the  results  of  the  30  runs.  Each  line  of  the  table  shows  the  values  of  the 
three  MOEs  (AOU  size,  mean  missed  distance,  and  accuracy)  that  were  computed  for  Nodestar  and 
MTST  for  the  indicated  run.  Each  run  covered  a  period  of  20  hours.  Every  15  minutes,  each  tracker 
produced  an  estimate  of  the  target's  location  for  a  total  of  81  times  at  which  estimates  were  produced 
(including  the  initial  estimate  at  time  0).  The  better  of  the  two  numbers  in  each  category  for  a  run  is 
labeled  with  an  asterisk.  A  quick  glance  shows  that  Nodestar  wins  most  of  the  comparisons.  The 
columns  labeled  G  and  H  give  the  number  of  times  at  which  contact  was  gained  and  held  respectively 
during  a  run.  The  sum  of  the  gains  and  holds  yields  the  total  number  of  contacts  made  by  all  three 
arrays  during  the  run. 

The  first  row  of  Table  B4  shows  the  results  for  the  case  of  a  convergence  zone  environment, 
constant-course-and-speed  target,  and  correct  estimates  by  the  Nodestar  of  the  FOM  and  propagation 
loss  curve.  Note  that  Nodestar  has  a  lower  average  AOU  size,  smaller  MMD,  and  higher  accuracy 
than  MTST.  The  next  two  lines  show  the  results  for  two  cases  with  the  same  target  motion  and 
environmental  characteristics  as  the  first,  but  with  Nodestar  assuming  that  the  FOM  is  one  standard 
deviation  higher  and  lower  respectively  than  it  actually  is  in  the  simulation.  Even  with  this  mismatch 
in  acoustic  assumptions.  Nodestar  wins  in  all  three  categories. 

One  would  expect  Nodestar  to  outperform  a  Kalman  filter  such  as  MTST  in  a  convergence  zone 
environment.  The  results  in  Table  B4  bear  that  out.  In  27  comparisons,  there  are  only  three 
instances  in  which  MTST  has  a  better  MOE  than  Nodestar.  Overall,  Nodestar  is  clearly  superior. 

In  the  case  of  a  direct  path  or  ducting  environment,  one  might  expect  that  Nodestar's  edge  over 
MTST  would  diminish  because  the  acoustic  information  in  this  environment  would  not  be  as  valuable 
as  in  a  convergence  zone  environment.  Surprisingly,  the  results  in  Table  B4  show  that  Nodestar 
continues  to  outperform  MTST  in  direct  path  and  surface  duct  environments  even  in  the  face  of  mis¬ 
estimates  of  FOM.  Examining  the  three  runs  (DCI,  DEI,  and  DPI)  in  which  the  propagation  loss 
curve  was  convergence  zone  but  Nodestar  thought  it  was  direct  path,  we  see  that  Nodestar  bettered 
MTST  in  eight  of  the  nine  MOE  comparisons. 
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Table  B4  —  Summary  Results 


Run 

AOU  (nm2) 
Nodestar  MTST 

MMD  (nm) 
Nodestar  MTST 

Accuracy  (%) 
Nodestar  MTST 

G 

H 

ccc 

2,839* 

4,880 

36.83* 

64.76 

69.76* 

27.62 

4 

13 

CCH 

2,354* 

4,158 

33.47* 

54.91 

79.91* 

33.36 

5 

16 

CCL 

946* 

1,122 

23.41* 

28.74 

70.19* 

45.18 

10 

49 

CEC 

1,571* 

2,530 

25.62* 

33.39 

67.19* 

65.17 

6 

52 

CEH 

2,874* 

5,822 

41.61* 

69.11 

52.34 

53.88* 

7 

8 

CEL 

2,908* 

4,641 

46.84* 

46.94 

46.38 

59.28* 

4 

23 

CPC 

4,597* 

7,584 

41.07* 

45.31 

67.58 

90.45* 

3 

9 

CPH 

4,017* 

7,247 

32.59* 

49.24 

88.3* 

65.90 

3 

19 

CPU 

11,080* 

16,860 

56.83* 

65.19 

84.48* 

81.9 

2 

0 

DCC 

1,057* 

1,743 

21.73* 

48.21 

50.7* 

30.39 

5 

32 

DCH 

4,852* 

8,849 

66.21* 

74.38 

19.04 

41.13* 

2 

7 

DCL 

1,911* 

3,620 

36.07* 

58.60 

25.61 

27.99* 

3 

26 

DCI 

1,482* 

2,595 

21.45* 

50.62 

63.57* 

36.42 

6 

25 

DEC 

3,705* 

6,688 

40.86* 

63.76 

71.65* 

61.95 

3 

5 

DEH 

2,012* 

3,697 

24.07* 

32.63 

71.41 

79.19* 

4 

29 

DEL 

3,768* 

6,147 

66.64* 

96.15 

30.05* 

25.37 

2 

13 

DEI 

4,664* 

8,129 

66.68* 

70.99 

27.02 

44.83* 

4 

5 

DPC 

898* 

1,322 

21.55* 

28.90 

70.83* 

68.99 

9 

24 

DPH 

2,155* 

3,831 

36.98* 

41.28 

41.42 

70.93* 

6 

16 

DPI 

2,670* 

4,706 

26.82* 

39.80 

72.15* 

61.60 

4 

26 

DPL 

6,680* 

10,940 

53.20* 

60.94 

62.53* 

58.07 

2 

0 

see 

5,668* 

8,162 

72.26 

63.92* 

39.96 

54.45* 

7 

3 

SCH 

1,224* 

1,529 

22.34* 

39.37 

83.21* 

40.36 

8 

30 

SCL 

5,755* 

9,273 

70.14* 

75.94 

53.64* 

43.45 

3 

4 

SEC 

24,920 

21,800* 

181.50 

133.10* 

9.77 

32.00* 

0 

0 

SEH 

3,275* 

5,139 

30.06* 

65.81 

89.13* 

50.01 

4 

22 

SEL 

3,160* 

4,410 

41.55* 

45.08 

65.58* 

65.25 

5 

16 

SPC 

1,750* 

2,587 

27.43* 

34.01 

71.12* 

70.91 

8 

13 

SPH 

1,639* 

2,310 

29.66* 

40.33 

62.48* 

50.30 

6 

20 

SPL 

860* 

1,261 

15.97* 

33.90 

89.47* 

43.12 

11 

22 

*  indicates  the  winner. 

Figures  B4  through  B6  graphically  present  the  comparisons  in  the  table.  The  comparison  of 
AOU  sizes  is  shown  in  Fig.  B4.  The  cases  are  labeled  along  the  horizontal  axis.  For  each  case,  a  bar 
shows  the  difference  between  the  AOU  size  for  Nodestar  and  MTST.  The  height  of  the  bar  indicates 
the  magnitude  of  the  difference.  A  black  bar  indicates  that  Nodestar  had  a  smaller  AOU  size  than 
MTST.  A  white  one  indicates  that  MTST  had  a  smaller  AOU  size.  One  can  see  that  Nodestar  beat 
MTST  in  all  cases  except  one. 
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Case 


Fig.  B4  —  AOU  comparison 

Figure  B5  shows  the  mean  missed  distance  comparison.  Nodestar  beats  MTST  in  all  cases  except 
two.  Figure  B6  shows  the  accuracy  comparison.  Nodestar  wins  21  out  of  30  of  these  comparisons. 

B5.  CONCLUSIONS 

The  obvious  conclusion  is  that  Nodestar  performs  much  better  than  MTST  in  favorable  situations 
(e.g.,  convergence  zone  environments  with  correct  acoustic  estimates)  and  continues  to  outperform 
MTST  in  cases  where  Nodestar  is  given  a  poor  estimate  of  FOM  or  even  the  wrong  propagation  loss 
curve.  The  conclusion  was  tested  statistically  using  a  Wilcoxen  test  for  paired  observations.  This  test 
showed  that  Nodestar  performed  significantly  better  than  MTST  for  all  three  MOEs. 

One  might  ask  whether  the  difference  between  Nodestar  and  MTST  could  be  attributed  to  the 
possibility  that  MTST  is  not  a  good  implementation  of  a  Kalman  filter.  In  our  opinion,  the  answer  to 
this  question  is  clearly  no.  MTST  is  one  of  the  best  Kalman  filters  for  submarine  tracking.  It  has  a 
motion  model  that  provides  sensible  velocity  distributions  and  it  handles  linear  information  in  a 
reasonable  fashion  given  the  linear  constraints  imposed  by  the  Kalman  methodology. 


Difference  in  Mean  Missed  Distance  (Nodestar  -  MTST) 


Technical  Documentation  of  Nodestar 
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Case 

Fig.  B6  — Accuracy  comparison 
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When  the  output  from  Nodestar  and  MTST  are  viewed  side  by  side  at  each  time  period,  the 
reasons  that  MTST  performs  so  poorly  compared  to  Nodestar  become  clear.  First,  the  solution  for 
range  information  in  the  bearing  detections  used  by  MTST  makes  it  difficult  for  MTST  to  localize 
the  target.  The  linear  approximations  required  to  incorporate  line-of-bearing  detections  cause  the 
target  distribution  to  get  drawn  into  the  sensor  in  a  sort  of  fatal  attraction.  Inability  to  use  negative 
information  allows  the  target  distribution  to  expand  into  regions  over  long  periods  when  no  contact  is 
obtained.  All  of  these  problems  will  be  shared  by  any  Kalman  filter.  They  are  not  specific  to  MTST. 
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